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INTRODUCTION 

The research which has been conducted during the past year under 
NASA Contract NAS8-30507 has consisted of two separate investigations, 
which, while related, are relatively independent of one another. The 
first of these consisted of a theoretical study of the tracking errors 
induced in an optical radar by atmospheric turbulence. Expressions for 
the magnitude of atmospherically induced angle of arrival fluctuations 
have been derived from Tatarski's theory of the propagation of light 
through randomly inhomogeneous media. These expressions have been com- 
pared with the available experimental data and found to agree within 
the experimental accuracy. We have also proposed a tracking system 
using multiple receivers which would reduce the apparent angle of 
arrival fluctuations without requiring large receiving apertures. The 
performance of such a receiver has been investigated theoretically 
and predictions of the RMS angular errors due to the earth’s atmosphere 
have been made. 

The second investigation consists of the reduction and analysis 
of certain government-supplied data. This data concerns the scintilla- 
tion of a 10.6 micron laser beam and the signal- to -noise ratio in a 
frequency modulated laser communications system. Techniques for 
reducing the subject data have been developed and software generated 
to implement the analysis on the IBM system 360-50 computer. The re- 
sults of the analysis are not complete at this time, however. Pre- 
liminary results will be presented. 
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ATMOSPHERICALLY INDUCED ERRORS IN AN OPTICAL TRACKING SYSTEM 
Introduction 

We have investigated the errors introduced into a laser tracking • 
system by atmospheric turbulence. The system being investigated is 
identical to those being used by Marshall Space Flight Center personnel 
to investigate atmospheric turbulence and is similar to systems which 
are proposed for early launch tracking of the Saturn V vehicle and for 
pointing deep-space optical communication systems. Theoretical ex- 
pressions for the magnitude of the root mean square tracking errors 
(angle of arrival fluctuation) as a function of the aperture of the 
receiving optics have been derived. Two assumptions have been made 
in deriving these expressions, viz; (1) that the tracking system effec- 
tively measures the centroid of illuminance in the focal plane of a 
well-corrected lens , and (2) that the fluctuations in the amplitude of 
the incoming wave may be neglected compared to fluctuations in its 
phase. 

An optical tracking system which was developed for Marshall Space 
Flight Center by Sylvania Electronics Corporation under Contract 
NAS8-20841 is currently being used by Marshall Space Flight Center 
personnel to investigate the effects of atmospheric turbulence. This 
system has been fully described in the literature [1], Without be- 
coming involved in the details of the operation of this tracker we may 
state that systems consist of receiving optics which focus an incoming 
laser beam onto the face of an image dissector tube. The image dissector 
measures the position of the focused spot as it moves across the face 
of the tube due to atmospheric turbulence. The output signal from the 
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system is proportional to the x and y coordinates of this spot. This 
data is then analyzed to yield, among other things, the root mean square 
deviation of the apparent angle of arrival of the laser beam in terms 
of the azimuth and elevation angles. We shall derive a theoretical 
expression for the RMS deviation in the azimuthal angle. A similar 
analysis holds for the elevation. 

The problem of defining the angle of arrival of distorted elec- 
tromagnetic waves is discussed by Beckman [2], Several authors have 
investigated the angle of arrival fluctuations of a laser beam propagated 
through the atmosphere using different definitions of the angle of 
arrival. Fried [3] has defined it as the angle of inclination of the 
plane which best fits the incoming wave front in a least mean square 
sense, while Heidbreder [4,5] has considered the problem in terms of 
the direction of maximum power radiation from a circular aperture 
antenna in a turbulent medium. Clearly for practical purposes one must 
take as the definition of the angle of arrival the quantity which is 
actually measured by the receiving system. In order to determine this 
quantity a careful analysis of the operation of the tracking system is 
required. If, for example, the system in question employs a quadrant 
detector with the output taken as the difference in signal from 
opposite quadrants then the apparent angle of arrival of the incoming 
beam as measured by the tracker will be proporational to the difference 
in illuminance on the two quadrants of the detector. For the case of 
a tracking system using an image dissector tube the situation is 
somewhat more complicated since the output will depend to a large ex- 
tent on how the signal is processed. We will assume that the tracker 
measures the instantaneous position of the centroid of illuminance in 
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the plane of the detector. While this assumption has not been fully 
verified by a detailed analysis of the particular tracking system in 
question it is felt that it should be sufficiently accurate for our 
purposes . 


Centroid of Illuminance 

Let x and y be Cartesian coordinates in the plane of the detector, 

x being taken in the horizontal direction and y vertically. With the 

assumptions discussed above the instantaneous azimuthal tracking error 

may clearly be expressed in terms of the moments of the intensity dis- 

\ _ 

tribution on the face of the image dissector tube by 9 x (t) = x/f where 
x is the x-coordinate of the centroid of the illuminance I(x,y) and f 
is the focal length of the system. 


_ 1. //xl (x,y)dydx 
f // I(x,y) dydx 


( 1 ) 


We have assumed that in the absence of atmospheric disturbance 
the centroid of illumination is at the origin. The limits of integra- 
tion will be assumed to be over the entire plane unless otherwise 
stated. I(x,y) may be expressed as the absolute square of the complex 
electric field quantity u(x,y). 


I(x,y) = u(x,y)u*(x,y) 


( 2 ) 


For a well corrected optical system u(x,y) is given in terms of the 
electric field in the entrance pupil U(£,q) by the well known relation 
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u(x,y) = A 


2xi 


U(s,n) exp { (x?+yn)} d^dn 


(3) 


A quadratic phase factor has been suppressed in equation (3) as it will 
not effect the illuminance. The integration extends over the aperture 
of the system. We may define U(£,q) to be zero everywhere except in 
the aperture and extend the limits of integration over the entire £-ri 
plane. 

We now let V = x/Af and W = y/Xf and combine equations (1) and (2). 


g ( t ) = x -^ v l u ( v > w )l dvdw 
//|u(v,w)[ 2 dvdw 


(4) 


From equation (3) we have 

-IJ 


j u(v,w) I dvdw 


*f 

x A 


IA f fu(^,n)e“ 2lri ^ v? ' H,n hcdn 
J . 

T1 i\ +2iTi(v? , +wr) 1 ) , , j 

U(? ,n )e d^'dn 3 dvdw 


(5) 


which may be rewritten 


j | u(v,w) | 2 dvdw = |a| 2 | 


f 


■ 

JJ 


J . 


uCc.n^a'.n') 


x exp {-2iti [v(?-£ ' ) + wCn-n^DdvdwdndCdn’d^ 1 


( 6 ) 


Performing the integration on v and w we obtain 


2 

u(v,w) | dvdw = | A 


rr 


U (£,q)U* (C 1 ,V) 


x 6 ')6 (n-n’ )d^dnd^ , dn' 


(7) 
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so that the integration on £' and n ' become trivial 


u(v,w) j 2 dvdw = | A| 2 JtJ (£,ri)TJ*(£ 1 ,n ' )d£dri 


Likewise the other integral in equation (3) is 


v|u(v,w) | 2 dvdw = |A| 2 Jj j U(?,n)U*(C' ,n') 

x v exp{-2iri [v(£-£ ' ) + w(n-*h , )3 dvdwd£dndc' dn* (9) 


On interchanging the order of integration this becomes 


f f C f 

v|u(v,w) | 2 dvdw = ] a| 2 U(?,n)U*(c' ,n') 

J J J 


x [ exp{-2Tfiw(ri-n T ) }dw] [ exp{-2Ttiv(^-C f ) }vdv]d?dc T dndri 1 (10) 


The integration on w leads to a delta function in the usual way. To 
perform the integration on v we note that 


6(5-5') - Je- 2 ’ i(S - ? ' )v dv 


Differentiating with respect to £ we have 


, „ 2 ,i e -2,iv(5-5') vdv 
3? 


2 7T 3 C 
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Mating use of this result, equation (10) becomes 


v|u(v,w) | 2 dvdw = | A | 2 | 


* 


. 

. 


u(?,n)U*a , ,n') 


X 6(q-n 1 ) 3lS f:i ^ ^d ridnMgdg’ 

ot, 


( 14 ) 


The integration of n is easily performed and the integration on t, may 
be carried out by parts. 



u(v,w) 


2 


dvdw = 


i 

2tt 



D(?,r^)U*(? , ,n , ) 


9? 


-d£d£ ' dri 


We take 


u = U(?,n') 


du 


3U 

9C 


dv . dc 

9? 


v = 6 (?-?») 


(15) 


(16) 


then 


U(c.n') = U(?,n’ ) 


f9U(t,n’) 

j 9e 


6(C“C')dC 


P(g,n , ) a5 ^~ g, '' t ig = U(c,n') 


’2 _ 9U(c' ,n 1 ) 

9? 


(17) 


(18) 


here £ a ^d are the values of £ at the edge of the aperture. In 
general for any aperture other than a rectangular one ^ and £ 2 will be 
functions of n. Substituting equation (18) into equation (15) we have 
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jv|u(v,w) | 


2^ j i A 
dvdw = 


2tt J 


u*(5’,n') 


x 


- uC^.n^fiC^-C 1 ) - 8U( 3g - ? T )t) d?’dn 


(19) 


hence 



u(v,w) dvdw 


2Trl A l 2 f i-l u <C 2 ,n) I 2 - ||0(Ei,n)l 2 
U*a,n) 3U 3 ( - ^ n) de dn 


( 20 ) 


The factor 1/2 in each of the first two terms arises from the fact that 
U(£,n) is discontinuous at and since it has been defined to be 
zero for all £ and n outside the aperture. Since integrals of the form 
/f(x)6(x)dx are properly defined only if f(x) is continuous at the 
point where the argument of the delta function vanishes we must exercise 
care in evaluating the integrals. First we perform a transformation of 
coordinates by letting £=£— 5 -^so that the argument of the delta function 
is zero at the origin. The limits of integration then extend from 
zero to infinity. We may now define U(-?) = U(£) rather than zero. 

This does not effect the value of the integral since U(-— £) is outside 
the range of integration. This integral may then be written as one-half 
the integral from minus infinity to plus infinity since the integrand 
is symmetrical. Hence the origin of the factor of one-half. 

Now we write U as luje 1 ^ so that the last term in equation (20 
becomes 


He 


-i* 


9 

9? 


U 


+i<j> j 
e v d£ 


( 21 ) 
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which, may be integrated by parts to yield 






2 dc + |iu(? 2 ,n) | 2 - f|ua r Ti) I 2 


( 22 ) . 


Substituting into (20) we have 


2 

jv|u(v,w) I 2 dvdw = “ j |u(£,n) | 2 d^dn (23) 


so that equation (4) becomes 


e x (t) 


A // i; |U( ? ,n )| 2 ||dcdq 

71 // z |u(c,n) I 2 a?dn 


(24) 


Equation (24) is a quite general expression for the instantaneous 
tracking error. By squaring and taking an ensemble average over all 
physical realizations of the incoming wave one would obtain the RMS 
angular fluctuation of 0 . Unfortunately, the resulting expressions 

X 

cannot be readily evaluated in terms of the statistical functions which 
are normally used to describe the fluctuations of the atmosphere. In 
one important case, viz, if the amplitude fluctuation can be neglected 
in comparison to variations in the phase, then equation (24) may be 
evaluated. We shall first turn our attention to this case and later 
discuss methods of evaluating equation (24) when this approximation is 
not valid. 


The Ray Optics Approximation 

Let us assume that the random fluctuations in the amplitude of the 
field may be neglected in comparison to the fluctuation in the phase. 
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This is essentially the ray optics approximation. The validity of this 
approximation has been widely discussed in the literature [6-8]. With 
this approximation ju(£,r|)| is a constant and may be taken outside of 
the integral sign in equation (24) . We then have 


e (t) 

x 


2tt // d£dn 


(25) 


where the integration extends over the aperture of the system. The in- 
tegral in the denominator is just the area of the aperture, hence 


e (t) 

X 


A 

2 it A 



C 1 

34 



d£dq 


(26) 


6 x (t) 


A 

2trA 


Pi 

[^> (^ 1> n) - 4>(4 2 » T i)]dn 

n 2 


(27) 


where and are th e appropriate limits of integration. We 

shall consider three cases: (a) a narrow slit of length a; (b) a 

rectangular aperture of length a and width b; and (c) a circular 
aperture of diameter D. 

(a) Narrow slit aperture: For this case equation (27) reduces 

to 


6 x ( t ) 


A 

2ira 


[*(■§) - <K- f ) ] 


( 28 ) 


Squaring 0 x 


and taking the ensemble average we obtain 






2 
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( 29 ) 


We recognize the quantity on the right as the phase structure function; 
hence the root mean square tracking error is 


6 rms 


A 

2ira 


^ D (a) ' 

<f> 


(30) 


Following Tatarski [9] we take D.(a) for a plane wave as 

<t> 


, _ 5/3 

D (r) » iTOS" ]|~| (31) 


for values of r much larger than the inner scale of turbulence £ . 

o 

Then 


o - A 'or _ -5/6 -1/6 

Q-q-m-o q » 0*00 IT 3. 

RMS 2tt o 


(32) 


For a<<£ Q , is proportional to r [10]. Clearly 6-^g must be in- 
dependent of the length for very short slits. From Tatarski' s 
equation 7.101 we have after some manipulation 


D = 4.21 1 _1/3 r ~ 5! 3 a 2 
<j> o o 


(33) 


Which yields a value 


9 SMS 


= . 326Ar 


-5/6 


o 


£ -!/6 


(34) 
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From equation (32) we find the value of (a) which will give this limit 
to be 5.2£ q . Hence we may write (32) and (34) as 


,429Jr ~ 5/6 a~ 1/6 

o 

a > 5£ 

o 

(35a) 

326Xr " 5/6 £ “ 1 ^ 6 
o o 

a < 5£ 

o 

(35b) 


(b) A rectangular aperture: Now consider a rectangular aperture 

of length a and height b. If we take the origin of our coordinate 

system at the geometrical center of the aperture, then = ~b, 
111 

^2 = ~ ’ ^1 = ~2 a anc * ri 2 = ~2 a * (27) then becomes 


e x (t) 


2irab 


t<Kc r n) 


<K? 2 >n)3dn 


(36) 


Squaring 6 and writing the square of the integral as a double integral 

X 

over ri and p ' we obtain 


6 


2 

x 



[<KC 1 ,ri) - <K£ 2 ,ri) 3 [$(£-£, h' ) 

~ <j>(4 2 ,n , )]dndp' 


(37) 


or 



x 



( 2ir ab J 




[<M4 1 >n)<j>(S 1 »n’) + (C 2 >n) (^ 2 ’ n T ) 


( 38 ) 


-4 | (4 1 J r i)<}>(? 2 ,P ) - (j)(c; 2 ,n)<{>(? 1 » r l , )]dndn , 
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2 

We now take the ensemble average of 0^ in equation (38) . Inter- 
changing the linear operations of averaging and integration we note 
that the resulting quantities in the right hand side of equation (38) 
are of the form as phase covariance functions. 

V*1 " ? 2> " < 39 > 

The phase covariance does not properly exist for nonstationary 
systems. We will introduce it, however, for mathematical convenience 
in manipulation without regard to the convergence of the integrals which 
are implied. will later be replaced by the structure function 
which is properly defined. The use of the covariance notation lends 
mathematical simplicity and leads to results equivalent to those which 
can be derived by manipulating equation (38) into the form of structure 
functions. Thus our procedure is justified even though it may lack 
certain mathematical rigor. Making use of the fact that (j> is a 
homogeneous and isotropic random variable, equation (38) can be written 
as 




’ X ' 

2 

n 2 

2irab 



J 

• 

n. 


[2C t (n-n'> 
$ 


1 "l 


(40) 


- 2C </ a^p^n-n'/mndn' 


By the well known relation between the covariance and structure func- 


tion, i.e.. 
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2C,Cr) = D M - D (r) 
4> <p 


( 41 ) 


equation (40) may be put in the form 


<o- 


X ] 

2 

z. 

[ n M 

2irab 



i 


■ 


hi 


Id (✓ (c 1 -c 2 r + (n-n’) 


T \2 


- D (n-n') 

9 


dqdr) ' 


(42) 


X 

2 -5/3 

6.88r ' 

o 

• n 2 

2irab 






'2 5/6 

[ (( ?1 -C 2 ) + (n -n’) z ) 


l ^1 


- I n-rT | 5/,3 3dndq’ 


(43) 


We first remove the absolute value signs from equation (43). The 
integrand in the first term is everywhere positive so that it does not 
present a problem. In the second term we may make a change of 
variables by letting x = q+a/2 and y = V+a/2. Noting that = a/2 
and n 2 ~ “ a /2, the second integral becomes 


r a 


r a 


1 5/3 , , 
x-y | dydx 


(44) 


Since the integrand is symmetrical about the line x = y and is positive 
for all x > y, we may write it in the following form which may be 
directly integrated. 


2 



(x-y) 


5/3 


9 11/3 

44 a 


dydx = 


(45) 
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The first term in equation (43) may be evaluated by expanding the 
integrand in a Taylor's series and integrating term by term. Noting 
that ?2 - ^l = we have 


T 2 2 5/3 

[(^-^r+Oi-rT) 3 dndn* = b' 5/J 


a/2 


-a/ 2 


_ t 9 5 / 6 

(l+C 3 -^-) ) dndn' 


+a/2 „ 


= b 


5/3 




. 


-a/2 l 


— (2 c , 4 c*7 ,6 

, , 5,n~n\ 5 ,n-n\ , 5*7,n-n \ 

1 + ) 9 — (t — ) + ~ 5 — (i; — ) 

6 b 6?2! b 6 ? 3! b 


5.7.13,n-n\ 8 . 

,4, , <■ b > + ' ‘ * 

6 *4! 


dndn ' 


(46) 


From equations (44), (46) and (43), we obtain 



S~6788 


X . 2 


(o=-> 


5/3 2ir 


- 2,-2 
a b 


9 11/3 . 2, 5/3 

-rr a + a b 
44 


+ 10 I (-i) n 


n 


2n+2 5/3-2n 

a b 1 


n=l 


6 n! (2n+l) (2n+2) 


(47) 


Where C = 1; 7; 7 x 13; 7 x 13 x 19; 7 x 13 x 19 x 25; etc. for 
n 

n = 1, 2, 3, ... . Equation (47) is valid for a less than b since for 
a greater than b the expansion in equation (46) does not converge. Now 
if we consider a rectangular aperture of length b in the direction 
along which the component of the tracking error is being measured and 
width a perpendicular to this direction, we define a shape parameter a 
for the aperture as the ratio of a to b. Clearly o must lie between 
zero and one for the series expansion to converge. Writing equation 
(47) in terms of a and taking the square root we obtain for 0 g 
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6 

rms 


4- /Xar r - 5/6 b" 1/6 
2lT o 


1 9 „ 5 / 3 

1 -“44 0 

2n 


+ 10 £ (- 1 ) 


C a 
n n 


1/2 


6 n n! (2n+l) (2n+2) 


(48) 


We see that for a rectangular aperture the RMS tracking error varies 
with the b dimension (i.e., the direction along which the component of 
angle is measured) like the reciprocal one-sixth power. The dependence 
on the other dimension is very slight since the term in braces varies 
from unity for a equal to zero (a narrow slit) to 0.96405 for a equal 
to one (a square aperture) . Clearly equation (48) agrees with our pre- 
vious results for the case of a narrow slit. 

We have evaluated the numerical coefficient in equation (48) for 
several intermediate values of the shape parameter. These results are 
given in Figure 1. 

(c) Circular Aperture: For a circular aperture of radius R, 

equation (38) becomes 


6 x - ( Sa>' 


R 

[ (C-L>n) d> ( c£,ti * ) + q>(c 2 » r i) ( l > (C 2 > 1 o r ) 

« 

-R 


- - 4>(C 2 » r >)^(C-^»n T ) ]dqdn (49) 


2 

Where A = irR is the area of the aperture and 


o 2 

q = (R -n ) 


1/2 


2 2 

? 2 = - (R -n ) 


1/2 


q = (R 2 -n 1 2 ) 


1/2 


q = - (R 2 -n’ 2 ) 


1/2 


(50) 



L .2 .3 .4 .5 .6 .7 .8 .9 1.0 

Numerical coefficient of equation (48) for a square aperture 
versus the shape parameter cr = a/b. 
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We now take the ensemble average of equation (49), making use of 
equations (31), (39), (41) and (50), and the homogeneity and isotropy 
of the structure function. We introduce new variables x = R and 

n 

y = ’ and obtain, after considerable manipulation 


< 0 - 


+1 +1 


2-nR 2) 


6.88 „ll/3 

5/3 R 
r 

o 


-1 -1 



/ 2 

1 2 

r i 

l' 

+ V 

t 1— / 

+ 

x-y 





>. j 



5/6 


£ 7') 2 




x-y 


5/6 


dxdy (51) 


The root mean square tracking error is easily obtained by taking the 
square root of equation (51). Introducing D = 2R, the aperture diameter, 
and simplifying somewhat, we obtain: 


6 


/ 6.88' X 


rms 5/6 2tr 
o 


2 7/12 _ 

— (J - J ) 

TT 


1/2 


D" 1/6 (52) 


where 


+1 +1 


'*■ | | t- 


xy± Jil-x') (1-y^) 


2 W , 2 , 


5/6 


dydx (53) 


-1 -1 


is a numerical coefficient which exactly corresponds to the shape co- 
efficient which was introduced for the square aperture. The integrals 
of equation (53) were programmed for the IBM 360 computer. Using a 
200 point Simpson's rule integration this yielded a value of .986 for 
the shape coefficient (the term in equation (52) contained in curly 
braces) . Thus we see that the tracking error for a circular aperture 
corresponds very' closely to that for a rectangular aperture with length 
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to width ratio of 1/2. 


Circular Aperture of Very Small Diameter 

For apertures which are small compared to the inner scale of 
turbulence the structure function is given by equation (33) instead of 
equation (31) . Thus for small apertures equation (51) becomes 



, A, 4.21 „4 

V 5 / 3 0 1/3 K 

IT 

o o 


+1 +1 


-1 -1 L 


[(/ 1-x 2 + / 1-y 2 ) + (x-y) 2 ] 


- [(/ 1-x 2 - / 1-y 2 ) + (x-y) 2 ] 


dydx (54) 


J 


which reduces to 


+1 +1 




4.21 


r 5 'h 1/3 
o o 


•16 


(1-xy) / (1-x 2 ) (1-y 2 )' dydx (55) 


-1 -1 

Since the x and y integrations may be separated we write equation (55) 


as 


<o- & 


67.36 


21T' r 5/3 £ 1/3 
o o 


r+1 


(/l-x 2 ) 


dx 


-1 



C 1 

/ — 


x > 

f 1-x 2 dx 


-1 



(56) 


The integrals in equation (56) may be evaluated directly yielding values 
of 2 and zero for the first and second terms respectively. Taking the 
square root to obtain the root mean square angular fluctuation we have 
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_ A / 67.36 
rms rr 5/ 6 1/6 

X" Af 


( 57 ) 


Thus we see that for a circular aperture the angular fluctuation 
will be independent of aperture size for apertures which are small 
compared to the inner scale of turbulence. 


Comparison With Calculations Based on Average Wave Front Tilt 

Fried [3] had developed expressions for the average tilt of a wave 
front after passing through a turbulent atmosphere, by expanding the 
phase <j> m a series of Modified Zemike Polynomials - i.e. , 


where 


Kx-jy) = a n F n (x l y) 


(58) 


q - C»rV 1/2 

F 4 = (irR 6 /12) (x 2 +y 2 -R 2 /2) 

F 2 = (ttR 4 /4)“ 1/2 x 

F 5 = (ttR 6 /6) (x 2 -y 2 ) 

F 3 = (TrR 4 /4)“ 1/2 y 

F 6 = ^ R6 / 24 )xy 


(59) 


Here a^ represents the average phase shift of the wave front, a^ and 
a^ its average tilt, a^ t ^ e spherical deformation, a^ and a^ the astig- 
matic deformation, etc. The average tilt is given in terms of 


a^ and a 2 by 


< a L/'< a l> +< ( 3 2^ 


(60) 


Fried has derived an expression for <^"a 2 ^ by requiring that the 
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polynomial expansion of <f> (equation (58)) gives a best fit to the 

actual phase in a least mean square sense. We have found that Fried's 

2 -1 

work contains an error in that he has omitted a factor of (ttR ) in 
his equation 4.6a. Following through Fried's work we find that the 
correct expression for a^^> is 



2 D 

= .883ttR (— ) 
r 

o 


5/3 


(61) 


The reader should compare this expression with equation 7.8a in Fried’s 
paper. 

With this correction in hand, we may proceed to calculate the apparent 
angle of arrival of the beam. Consider a plane wave traveling in approxi- 
mately the z direction. Neglecting deformation of the wave front we 
may write its phase as 




a 2 F 2 + a 3 F 3 


(62) 


2irz 

But the phase is also (— — ) so that the equation of the isophase 

A 

surface is 


2ttz 

X 


+ a 2 F 2^ + a^F^(y) = constant 


(63) 


We let K be a vector normal to the isophase surface and k. be the unit 
vector in the z direction. Substituting the values of F from equation 
(59) into equation (63) we may write K as 
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K = GRAD 


z + 


,A 2 


(a 2 x + a^y) 


(6A) 


Since 


|K| cos0 t = K . K 


(65) 


a straightforward calculation yields 


cos0 T 


2 2 l~l/2 

TT (a 2 + a 3 > + 


hr R 


J 


( 66 ) 


Squaring both sides and noting that for small angles 


COS0 T = 



(67) 


and 


i . * / 2 . 2. 

1 + 34 (a 2 + a 3 5 

ir R 


r l/ 2 


= 1 v(a 2 + a„ 2 ) + ... (68) 

2ir R 


We obtain 


SMS 




170 J <Ay 

IT R 


(69) 


Substituting from equation (61) 


RMS 


A- D - I/6 


5/6 


.414 


(70) 
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Chase [11] has noted that the constant (.883) in equation (61) 
should properly be .828 since there were numerical errors in Fried's 
claculations which were later corrected by Chase. 

Heidbreder’s calculations [4,5] based on the maximum of the radiation 
pattern of a circular aperture also led to a minus one-sixth power de- 
pendence on the aperture diameter but with a somewhat different proportion- 
ality constant. That the calculations of Fried, those of Heidbreder 
and ours all lead to the same functional dependence with only slightly 
different numerical constants is consistent with the conclusion (as 
Fried and others have pointed out) that wave front tilt is the pre- 
dominant effect responsible for angle of arrival fluctuations. 

Comparison With Experimental Results 

Wyman [12] has reported experimental measurements of the RMS track- 
ing errors as a function of receiver aperture over 3.2 and 6.4 Km paths. 
Figure 2 shows the root mean square angular fluctuations for both path 
lengths and several aperture sizes . Theoretical curves of angular 
fluctuations versus aperture size from equation (52) are also plotted 
in Figure 2. Here the correlation distance r^ has b,een chosen to pro- 
duce the best fit with the experimental data. 

The amount of data available is very minimal so that the theory 

cannot be said to be fully verified. Nevertheless the available data 
-1 j c. 

does fit a D power law to within the experimental error. 

TRACKING SYSTEM USING MULTIPLE APERTURES 


From our previous results and from the available experimental data 
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SMS Fluctuations (p.rad.) 



Aperature (cm) 


Figure 2. RMS angular fluctuations in microradians as a function of 

receiving aperture diameter for 3.2 and 6.4 kilometer paths. 
Experimental data after Wyman, Reference 12. 
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it is clear that when one increases the receiving aperture of an optical 
tracking system the atmospherically induced errors are decreased due 
to averaging over a larger part of the incoming wave front. Since the 
critical factor in reducing the atmospheric errors is the linear 
dimensions involved, it seems likely that an effective increase in 
accuracy could be obtained by using two tracking systems with small re- 
ceiving apertures located some distance apart instead of a larger single 
aperture system. Since large diameter optics are extremely expensive 
the two aperture systems might well be more economical than a larger 
single system. 

For this reason we have studied the error in a double receiver 
tracking system. In the process of our investigation we have devised 
an experiment to test the validity of our theoretical calculations which 
should overcome some of the difficulties previously encountered in 
measurements of atmospheric effects. This experiment will be discussed 
later. 

Theoretical 

Let us consider two identical ^optical tracking systems, such as have 
been previously described, located in close proximity to each other and 
observing a common target. Each tracker will see an apparent angular 
motion of the target due to atmospheric fluctuations which we will de- 
note by 0^(t) and ©^(t) respectively. Let us suppose that we may observe 
either the "relative tracking error", 0^(t), which we shall define as 
the instantaneous difference in 0^(t) and 0 2 (t), or the "average track- 
ing error", 0 (t) , which we shall define as the average of 0 (t) and 

3 - 

0 2 (t) . As shown in Figure 3 we take E, and p to be cartesian coordinates 
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in the plane of the receivers' apertures and let R be the radius of each 
aperture and L the distance between their centers. If 0^ and 0^ are the 
components of the instantaneous angular tracking errors of each system 
along the n axis then, as has been shown. 


- "2wA 


[$(5 1 ,n,t) - K£ 2 ,ri,t)3 dn 




e 2 (t) 


2irA 


[$(C 1} ri,t) 


$(S 2 ,n,t)] dn 


(71) 


(72) 


where n^, n 2 » n^, and n 2 are the appropriate limits of integration for 
a circular aperture, and £ are functions of n, and A is the area of 
a single aperture. 

We have defined the relative tracking error as 


6 = 0 - - 
r 1 2 


(73) 


so that 


0 = 


r 2irA 


t$(5 15 n,t) - $(E 2 ,n,t)] dn 


rh-i 


Tit 


[ (C 1 ,n.,t) - K? 2 ,n,t)] dn 


(74) 


Making use of the fact that the centers of the apertures are located 


at (0,± respectively, we obtain the relations 
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= / R 2 -(n - |) 2 


52 = -/R 2 -(n-|) 2 




- /„2 , , L .2 

«2 = - / R _(„ + j) 


n l 2 


"2=f +S 


2 - R "2- 


Squaring equation (74) one obtains 


9 >> ' 


f n 2 f n 2 


t) 


§(C 2 ,n,t) 


n l n l 


x *(€j_jn f ,t) - ,t)j dqdn 


rn 9 rn 9 _ _ 

+ tUj.n.t) - $(€ 2 »n,t) 


n l n l 


X $(T^,n’,t) - $(Ij,n',t) dridn’ 

’b 9 rb 2 

,n,t) - 4(5 9 ,n,t) 

J r\ r\ . 


"1 n i 


X $Uj_,Ti',t) - 'J(5 2 ,n , ,t)l dndn' 
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■ n 2 

rh 0 

2 

J 

hn 


<Kq,n,t) - $(q,n,t) 


x 


< Kq,n' ,t) -<Kq,n',t)] dndn’ (77) 


We expand the products under the integral sign in equation (77) and 
take the ensemble average of both sides. The 16 terms resulting on 
the right hand side are recognized as phase covariance functions [13] . 
Using the procedure we have derived in the first section we replace 
these by appropriate phase structure functions using the relation 


2C „ - V”> - V r > 


(78) 


or 


V r) - 2C t (0) - 2C 4 (r) 


(79) 


Performing these operations and substituting from equation (75) we 
obtain for the mean square of the relative tracking error 


<6 





q[/(q - V Z ) Z + (n-n ') 2 ' 



q ) 2 + (q-q ') 2 ‘ 
1 J 


dqdri T 
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(80) 


+D *K (5 1 " + ^-n') z j 

As before we take the structure function as 


-i a . i 


dridn' 


D (r) =6.88 

<P 


(81) 


where r Q is the wave correlation distance. Combining equations 75, 80 
and 81 and making a change of variables in each integral of the form 


i , . L. 

x * e (n± i> 


(82a) 


1 . L, 
y - 5 (n± ~) 


(82b) 


we have 
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Here we have taken, the signs in equations (82a) and (82b) as (+,+) 
respectively in the 1st integral, (-,-) in the second, (-,+) in the third 
and (+,-) in the last. 

Equation (83) may be put into the form 


<e 2 > 

r 


tO 2 6.88 R 11/3 

M r 5/3 




(0) - I + (L/R) + I (L/R) ] 
P P 


(84) 


where 


I“(L/R) 


r+1 r+1 

-1 J-l 



+ (x-y+L/R) 2 


5/6 


dxdy 


(85) 


Comparing this to our previous expression for the tracking error due to 

2 

a single circular aperture, <0 >, we see that the first two terms in 

2 

equation (84) are just 2<0 >, while the last two terms are the reduction 

X 
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in error obtained by using a two aperture system, i.e., 


<0 


2 

2<e > 

X 


fx l 2 6.88 R 11/3 

,+ (L/R) - r (L/R) 

| 2 ttaJ 5/3 

J p p 


( 86 ) 

1 


o *- 

Examination of equation (84) shows that for L equal to zero (the two 

2 

systems sharing a single aperture), c0 > is zero. This is clearly 

necessary since for this case (if it were physically realizable) 0^(t) 

4* _ 

is identical to 0„(t) for all t. For very large L, I (L/R) = I (L/R) , 

^ r F 

and equation (86) reduces to 


2 2 

<0 > = 2<0 > 
r x 


( 87 ) 


or 


0 (Double Aperture) = /T 0 (Single Aperture) (88) 
rms rms 


We next turn our attention to the component of angular tracking 
error in the 5 direction (i.e., at right angles to the line of centers 
of the two apertures). An analysis similar to that performed above leads 
to an equivalent expression for this component, i.e.. 


<0 > « 2 
r 


12 6.88 R 11/3 * 


2 it A 


5/3 


I + (0) - T(0) - I* (L/R) + l"(L/R)j 
s s s s 


(89) 


where I (L/R) is now given by 

r+1 r+1 -r 


i;(l/r) = 


-1 


-1 l 


---y ± / L _ y 2 + h 


+ (x-y)' 
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dydx (90) 
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As before this reduces to zero for L equal to zero and approaches the 
square root of twice the RMS error of a single aperture for large L. 

From equation (84) (or (89)) we may find the ratio of the relative 
tracking error between two apertures to the RMS error in a single 
aperture system of the same diameter, i.e., 

1/2 

(91) 

where I ± (L/R) denotes either I (L/R) or I (L/R) (given by equation (85) 

p s 

or equation (90)) depending on whether the angle being measured is along 
the line of center or at right angles to it. 

The quantities I - (L/R) have been evaluated on the IBM 360 computer 
using a 201 point Simpson's Rule integration for values of L/R from 2 
to 500. The results are given in Table I and Figure 4. 

We have defined the average tracking error 6 (t) as 

cl 

e n (t) + e (t) 

e a (t) = — 2 ~ (92) 

The RMS value of this quantity is easily found by a straightforward 
modification of the previous calculation. 

1/2 

(93) 

This quantity has also been evaluated as a function of (L/R) , the re- 
sults being given in Table I and plotted in Figure 5. 


A_ _ JL I (0) - I (0) + I (L/R) - I (L/R) 

"2 2 + - 

:> i (o) - 1 (o) 
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Discussion of Results 

A double tracking system such as has been described might operate 

by taking the average of the output of the two trackers as the actual 

position of the target. From Table I it can be seen that for a 

separation of 10 radii between the centers of the two apertures, the 

2 1/2 

SMS tracking error <0 > is approximately 86% of the error for a 
single aperture and for a separation of 100 radii the error is 78%. 

Since the RMS tracking error for a single aperture varies as the 
reciprocal one-sixth power of the diameter we see that for 10 radii 
separation the error will be the same as that in a single system whose 
aperture is 2.5 times larger. That is a double system with two 20 cm. 
diameter receivers located 1 meter apart will have no more error than 
a single receiver with a 50 cm. diameter objective. Likewise if the 
receivers were 10 meters apart the system should be subject to tracking 
errors approximately the same as an 88 cm. diameter single receiver. 
Whether or not the gain in accuracy in a double receiver system would 
justify the additional complexity will probably depend on the require- 
ments of the particular system in question. 

Experimental Verification of the Theory 

One of the major problems encountered in the experimental investi- 
gations of atmospheric effects is that the statistics of the atmosphere 
are non-stationary. Thus we are always attempting to gather statistical 
data on a system whose statistical properties are constantly changing. 
For example, if one attempts to measure the angular fluctuations of a 
beam propagating through the atmosphere as a function of receiver 
aperture size he finds it difficult to collect enough data to be 
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TABLE I 


L/R 

AVERAGE ANGULAR 
FLUCTUATION 

RELATIVE ANGULAR 
FLUCTUATION 


PARALLEL 

PERPENDICULAR 

PARALLEL 

PERPENDICULAR 

2 

0.9553 

0.8980 

0.5912 

0.8799 

3 

0.9316 

0.8701 

0.7268 

0.9858 

4 

0.9147 

0.8546 

0.8084 

1.0385 

5 

0.9019 

0.8441 

0.8637 

1.0723 

6 

0.8919 

0.8362 

0.9045 

1.0967 

7 

0.8839 

0.8300 

0.9352 

1.1155 

8 

0.8769 

0.8249 

0.9613 

1.1305 

9 

0.8710 

0.8206 

0.9824 

1.1430 

10 

0.8659 

0.8169 

1.0003 

1.1536 

50 

0.8039 


1.1894 


100 

0.7849 


1.2392 


500 

0.7534 

0.7384 

1.3151 

1.3487 

999 

0.7441 


1.3361 



statistically significant in a short enough time to insure that the 
atmospheric donditions have not changed during the course of the experi- 
ment. The alternate approach of collecting data over a long period of 
time, say many months, and then assuming that the average results represent 
some sort of hypothetical "average atmosphere" is not only time consuming 
and expensive but also is unsatisfying since the variation in atmospheric 
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Figure 4. Ratio of Relate 
to Fluctuatior 





Figure 5. Ratio of Average Angular Fluctuations to Fluctuations in A Single Aperture 
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turbulence from day to day may be as large as or larger than the effect 
which one wishes to observe. Furthermore one is also tempted to discard 
data which does not agree with the theory on the grounds that "the 
atmosphere must have changed during the run." The validity of data sub- 
ject to such subjective interpretation must always be suspected. One 
would like an experiment which would confirm the essential validity of 
the theoretical model of the atmosphere yet not be dependent upon the 
strength of the turbulence at a particular time. 

We believe that a two aperture tracking system would provide such 
an experimental opportunity. The scale of turbulence and the turbulence 
strength enter into our calculations only through the single parameter 
r , the correlation distance. This parameter has been eliminated from 
equations (91) and (93). Thus an experiment in which the RMS relative 
and/or average angular fluctuations were measured simultaneously with 
the individual angular fluctuations in each tracker should provide a 
sensitive test of the theory independent of the value of r^. Any in- 
accuracy in the theoretical model of the atmosphere, such as a departure 
of the atmospheric fluctuations from a log-normal distribution, should 
be easily demonstrated. 

GAUSSIAN DISTRIBUTION OF BEAM INTENSITY 

We have also considered the effects of a Gaussian distribution of 
intensity of the laser beam on the apparent angular fluctuation. A 
distribution of intensity across the beam of the form 



o 


( 94 ) 
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was assumed. From equation (94) we have that 


ff |u(£.n ) | 2 d£ dn 

a on - x 1 3 ? 

e x (t) ~ 2m 


// |u(5.n) | d£ dq 


( 95 ) 


which yields 


„ e -°aV> || d5 d n 
x 

// E e 0(5 +n 5 d? dn 


(96) 


where a .is a parameter which describes the beam width and the inte- 
gration extends over the receiving aperture. The integral in the 
numerator of equation (96) may be integrated by parts on £ to eliminate 
the derivative 3$/3£. The resulting expression is then squared and the 
ensemble average taken in the usual way. After more manipulation the 
resulting expression can be recognized as phase structure functions. 
Unfortunately this expression contains 16 terms each of which is a 
multiple integral of the product of and an exponential function 
of the coordinates squared. Sufficient computer time has not been 
available to evaluate these quite complex expressions. 


EFFECTS OF ATMOSPHERIC TURBULENCE ON 
A C0 2 LASER COMMUNICATIONS SYSTEM 


Introduction 

The second phase of the research conducted during the course of this 
project involves the reduction and analysis of certain government 
supplied data concerning the effects of atmospheric turbulence on a CO^ 
heterodyne communications system. This data was collected during the 
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Summer of 1969 at Marshall Space Flight Center and was transmitted to 
the University of Alabama for analysis m September, 1969. We have 
generated the necessary computer software to reduce the data and have 
begun to analyze it. Due to the large amount of data and the limited 
time available this analysis is not yet complete. We expect, however, 
to be able to complete the analysis in the near future. In this re- 
port the problems encountered in analyzing the data are discussed; 
the computer programs which have been written are described; and 
preliminary results are given. 

Experimental 

The atmospheric experiments were performed by the project director 
in association with MSFC personnel at MSFC and are therefore not a part 
of this project. We will therefore not attempt a complete discussion 
of the experimental procedures in this report but will give only a brief 
description necessary to describe the data which was analyzed. A 
complete discussion of the experiment is being prepared by the project 
director and MSFC personnel and will be published in the near future 
as a NASA Technical Report. 

The subject data consist of measurements of the scintillation and 
heterodyne signal of a CO^ laser beam propagated over a 3.2 Km path 
between Madkin Mountain and the Astrionics Laboratory building, both 
located on Redstone Arsenal, Alabama. The laser transmitter was a 
2 watt, stabalized CO^ laser with a 10 cm., off-axis Cassegrainian 
collimator. The laser could be frequency modulated by driving one of 
the cavity mirrors mounted on a piezoelectric pusher. The receiver 
consisted of a 10 cm. aperture off-axis Cassegrainian telescope, a 
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local oscillator lasejr, combining optics, and a Hg-CdTe detector. The 
receiver was fitted with removable aperture stops so that its aperture 
could be varied from 2 cm. to 10 cm. The transmitter and receiver were 
constructed for MSFC by the Minneapolis-Honeywell Corporation and 
have been described in the literature [14] . 

Signal-to-noise measurements were made by extracting the 10 MHz 
beat note between the received signal and local oscillator after it had 
passed through one stage of IF amplifications. The signal was then 
detected with a simple diode circuit and the resulting voltage re- 
corded on magnetic tape. 

Scintillation measurements were made by turning off the local 
oscillator laser and chopping the transmitted beam at 90 Hz by means 
of a mechanical chopper located in the transmitter. The output of 
the detector was amplified and recorded directly on a 14 channel 
magnetic tape. Thus the recorded signal consisted of an amplitude- 
modulated square wave whose amplitude was proportional to the 
instantaneous power being received. In this way fluctuations in 
background illumination were eliminated. 

Twelve channels of each analog magnetic tape were used for data, 
the other two channels being reserved for identification and timing 
purposes. The tapes were digitized at a sampling rate of 1 KHz in 
order to reduce the time required for A/D conversion. The resulting 
digital tapes were recorded in a multiplex format with five channels 
of the analog tape on each digital tape. 

Between June 15 and August 31, 1969, approximately 750 observations 
of scintillation were made at all hours of the day and under as wide 



42 


a variety of weather conditions as possible. Each observation con- 
sisted of about 90 seconds of recorded data from which a 60 second 
segment near the middle of each run was selected for digitization. 

The time of day, temperature, humidity, wind speed and general 
weather conditions were recorded for each observation. 

Data Reduction 

In order to reduce the data a program has been written for the 
IBM 360 Model 50 computer. The principal problems which were en- 
countered in writing this program concerned formating the data for 
the computer and extracting the amplitude of the square wave. The 
latter proved to be somewhat difficult since the sampling rate during 
digitization could not be accurately synchronized with the period of 
the square wave. The sampling rate of 1 KHz and the chopping rate 
of 90 Hz should yield approximately 10 samples per cycle of the square 
wave. In actuality it was found that the number of samples per cycle 
varied between 10 and 12 due to the sampling rate not being an in- 
tegral multiple of the square wave frequency. It was therefore 
necessary to design a program which would determine whether a particu- 
lar data point was a base point (i.e., from the part of the square 
wave when the laser beam was blocked by the chopper) or a signal point 
(when power was being received from the laser beam) . The problem was 
further compounded by the fact that the rise and fall times of the 
square wave were non-negligible so that about one percent of the data 
points were sampled during the switching transient and should there- " 
fore be discarded. Furthermore it was found that some of the data 
contained an occasional noise spike which should be eliminated. It 
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was decided that the elimination of these. spikes would not adversely 
effect the validity of the analysis so provisions for eliminating them 
were also included in the program. 

The program which we have written to divide the data into base 
and signal points operates basically as follows. One record, contains 
ing 2000 characters, is read from magnetic tape. These 2000 characters 
represent 200 sample values from each of 5 data channels. Each sample 
value is a 10 bit binary number plus sign bit occupying two tape 
characters. The 400 characters corresponding to the channel being 
analyzed are converted to internal floating point format and are stored 
in an array. A second record is read from tape, converted, and stored 
in a second array. To begin the analysis twenty data points from the 
first of the array are selected and the maximum and minimum are found. 
Two limits, and L are then set by the relations 


L 

L 


1 

2 


A - P, (A - A . ) 
max 1 max mm 


A . + P 0 (A 

mm 2 max 


- A . ) 
mm 


(97) 

(98) 


where A and A . are the maximum and minimum of the first twenty 
max mm 

points and P^ and P^ are constants between zero and one-half. Since 
the signal was inverted when it was recorded on analog tape, the base 
line is greater than the signal, hence a particular point greater than 
is considered a base point, if it is less than 1^ it is considered 
a signal point. Points lying between and are assumed to be from 
the transient portion of the wave form and is neglected. 

The computer is programmed to take each point successively and 
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determine if xt is a base point, a signal point or neither. As a 
preliminary to processing the first twenty points are scanned and the 
begxnning of a base line segment of the wave form xs found. Then new 
limits are set on the next 15 points and they are scanned and grouped 
into 3 arrays, a base line segment, a signal segment and a second 
base line segment. Each array may contain up to 7 points. At this 
time the amplitude of the square wave is computed for the group of 
signal points (as will be described later) and stored. The second 
group of base points is transferred into the first array, new limits 
are set using the next 10 data points , a new group of signal points and 
base points are found to fill the second and third arrays, and finally 
their amplitudes are computed. This process is continued until the 
200 points from the first record have been used. At this time the 
200 points from the second tape record are transferred into the array 
formally occupied by the first record and a third record is read from 
tape. Processing then continues until all points in the data set have 
been processed. 

Data Processing Irregularities 

As previously mentioned, several irregularities in the data are 
possible and a number of checks have been built into the program to 
provide for them. These checks are as follows : 

1) During the search for either base points or 
signal points more than 10 consective points 
are found. 

2) After completing a search less than 3 base or 
signal points have been found. 

3) More than 3 consective points satisfying neither 
the base or signal point criteria are found. 
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Any of these three conditions indicate that the waveform is departing 
drastically from a modulated square wave and appropriate action should 
be taken. For conditions one and three the program skips 10 data 
points, or approximately one cycle of the square wave, and begins 
processing again. For condition two the computation of amplitudes 
are surpressed and processing continues. As a further check, i'f the 
total number of errors in any record exceeds five the entire record is 
omitted. 

During processing a record is kept of each time an irregularity 
was encountered and this information is printed in tabular form at the 
completion of processing. Clearly if an excessive number of irregulari- 
ties occurs in a given run the results of that run must be suspect. 

Computation of Amplitudes 

After each cycle of the waveform has been processed to divide the 
data points into base points and signal points the amplitude is 
computed. Three methods for computing the amplitude have been tried. 

The first method took the base line for a group of signal points as the 
average of all the points in the group of base points preceding it and 
the one following it. That is, the background during the half-cycle 
in which the laser beam intensity was recorded is taken as the average 
background recorded during the half-cycle immediately preceding the 
signal and the half cycle immediately following it. This average base 
line was subtracted from each signal point and the difference taken 
as the amplitude of the laser beam at that instance. 

The second type of amplitude calculation considered was to re- 
construct the base line during the period when the signal was recorded 
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by fitting a least-mean-square curve to the base line points on either 
side. This method was found to give very erratic results and was 
abandoned. 

The third method consisted of taking the difference in the first 
signal point in a group and the last base point preceding it as an 
amplitude. The difference in the last signal point in the group and 
the first base point following it gives a second amplitude. This 
method gives only two amplitudes per cycle but has the advantage that 
they are evenly spaced. 

The final computer program contained both the first and third type 
amplitude calculation, the one to be used being selected by a parameter 
read during execution. An option is also provided which will either 
store all amplitudes along with the time which the amplitude occurred 
or will, instead of storing the amplitude, count the number of times 
an amplitude lying in a given range occurs. The former yields re- 
ceived beam intensity as a function of time while the latter gives the 
probability distribution function for the intensity fluctuations. 

Either is saved for whatever analysis one wishes to perform on the data. 

Program Checkout and Adjustment of Parameters 

In order to test the program a segment of data was printed from the 
magnetic tape and was inspected. Each data point was classified as 
either a signal point, a base point, or a bad point from the transient 
portion of the waveform. This classification was purely subjective, 
yet in inspecting the data there was usually no questions as to how a 
particular point should be categorized. The same data was then fed 
into the computer. Print statements were added to the urogram to list 
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each, point and indicate how the computer classified it. The program 
was run several times varying the parameters P^ and P^ (equations (97) 
and (98)) between runs and the results compared with the subjective 
analysis. Data having "bad places" in it was also processed and the 
results compared with our judgement as to whether or not a segment 
should be omitted. On the basis of these comparisons the parameters 
P^ and P^ were set at 0.05 and 0.10, respectively. That is a point 
within 5% of the maximum base point or 10% of the minimum signal point 
would be retained while points between these limits were discarded. 

These limits seemed to allow the computer to make very nearly the 
same decisions as we would have made had we analyzed the data by hand. 

It is felt that due to the statistical nature of the analysis whether 
or not a few points are discarded as bad when they should have been 
retained will not appreciable effect the results. Preliminary analysis 
of the data seems to confirm that the results are not too sensitive 
to small changes in the values of the limits. 

Calculations of the Scintillation Statistics 

The final segment of the data analysis program accetps the proba- 
bility density function for the intensity fluctuations which has been 
generated and computes the scintillation statistics. The program 
computes and lists the class mark for the intensity and the corre- 
sponding value of the log-amplitude defined by 

= |(ln) I i /T (99) 

where and are the log-amplitude and the intensity for the ith 
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class interval, and I is the mean intensity. The frequency for each 
class and the cumulative probability are also listed. 

It has been customary in the literature to test the hypothesis of 
log-normality of scintillation data by plotting the cumulative proba- 
bility function of the log-amplitudes against a "probability scale" 
such that if the data is log-normal the resulting curve will be a 
straight line. Not only is this procedure not a very sensitive test 
of a statistical distribution but it is also very time consuming 
when a large quantity of data is to be processed. We have therefore 
included in the program a chi square test on both a normal and a log- 

t 

normal distribution function. These tests provide a quick and sensi- 
tive means of testing the hypothesis of log-normal scintillation. 

The statistical analysis routine also computes the mean, standard 
deviation, skewness, and kurtosis for both the intensity distribution 
and the log-amplitude distribution. From the standard deviation of the 
log amplitudes the atmospheric structure constant will be computed. 

The skewness and kurtosis are computed to give an additional check on 
lo g-normal r ty . 

Description of Program 

A listing of the computer programs described above along with a 
sample output are included in Appendix A, and a sample output is shown 
in Appendix B. For reference, Appendix C lists a program for generat- 
ing log-normal random numbers which was used in checking the chi 
square test subroutines. 


Atmospheric Structure Constant 
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The log amplitude variance C (0) for a plane wave is given in terms 
of the atmospheric structure constant by [15]. 

C„(0) = 0.309 k 7/6 Z 11//6 C 2 (100) 

where k is the wave number and Z the length of the path through the 
atmosphere. For a spherical wave the corresponding expression is [16] 

C„(0) = 0.124 k 7/6 Z 11/6 C 2 (101) 

t n 

In our preliminary analysis we have neglected the finite 

aperture of the receiver and treated it as a point source. Equations 

2 

(100) or (101) can then be used directly to obtain by noting that 
the standard deviation of the log-amplitude distribution which we have 
computed is the square root of C^(0). It is well known however that 
a finite receiving aperture has the effect of averaging the intensity 
fluctuation from various parts of the wave front thereby reducing the 
variance of the scintillation. This effect has been investigated by 
Fried [17]. From Figure 2 in Fried's paper we note that for large 
scintillation conditions and for the path lengths and apertures used 
in the experiment, this effect will be significant. 

In order to allow for aperture averaging we may use the expressions 
given by Fried [14,16], viz. 

2 

ff S 2= [f ° 2 0 C I (0) (102) 

2 

where is the signal variance which corresponds to the square of the 
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standard deviation of the intensity fluctuation, D is the diameter of 
the receiving aperture, and 0 is an aperture averaging factor given by 


0 = 


16 

•ri)" 


D exp [4C (p)]-l 

pdp exp I4C (0)]-l' H(P/D) 


(103) 


H(p/D) is the optical transfer function of a circular aperture 


H(p/D) = cos^Cp/D) - (p/D) [1 - (p/D) 2 ] 


1/2 


C104) 


and C (p) is the log-amplitude co-variance given by 


C,<P) = c ? (0) l 


n=0 


a + b 
n n 


kp_ 2 } 

4z J 


X 


kp 

4z 


/ (2n) 


- 7.53034 


kp 

4z 


5/6 


C105) 


In the last expression a^ and b^ are the expansion coefficients 
for the modified confluent hypergeometric function i x ) and 

are given by 


a = 1 b = 6.84209 

o o 

and the recursion relations 

a n = 'Vl [ (2n " 23/6 ^ 2n ~ 17/6)/(2n-l)(2n)| 

k n = -b n _ 1 |(2n - 17/6) (2n - 11/6) (2n-l)/(2n) (2n+l) 2 


(106) 

(107a) 

(107b) 
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With the additional relation that the intensity variance 0^(0) in 
equation (102) is related to the log variance by 


CjCO) = 


exp [40^(0) ] 


-1 


(108) 


2 2 

equations (102 - 108) specify C^(0) in terms of and I . Since 
I and a are iust the mean and standard deviation of the received 
intensity we may compute C (0) and then using equation (101) find the 
atmospheric structure constant. 

Combining the above equations we have 


r 


,2 


o 

s 




j^exp 




x 


fD 

^ o 


P dp 


exp [4f C £ (0)] -1 

exp [4 C £ (0)] -1 


H(pD) 


(109) 


where f(kp/4z) is the summation given in equation (105). 

Since a /I is an experimentally determined constant, equation (109) 
s o 

is an integral equation for C. (0) . We have programmed the IBM 360 computer 

to solve equation (109). The technique used is to evaluate the integral in 

equation (109) for a number of trial values of C^(0) using a fourth order 

Runga-Kutta integration. This gives a table of a / I as a function 

of C (0) . from this table the value of C (0) corresponding to the 
36 36 

measured value of a /I is determined using Lagrange-Hermite inter- 
polation formula. 


Preliminary Results 
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At the present time approximately 60 of the 750 scintillation 
measurements have been processed. This includes computation of the 
cumulative probability curve, the moments of the intensity dis- 
tribution, the moments of the log-amplitude distribution, and a 
chi square test for normal and log-normal distributions. Aperture 
averaging effects have not yet been considered nor has the power 
spectral density. 

From the data thus far analyzed several definite trends are 
apparent. It must be emphasized, however, that these results are 
based on a cursory inspection of a small part of the data so that 
any conclusions drawn at this time must be considered as only 
tentative. We shall, however, discuss some of these results briefly. 

1) Aperture Averaging . The scintillation data consist of groups 
of runs in which the receiver aperture was varied from 2 cm. to 
10 cm. The runs within a particular group were made in as rapid 
succession as possible in order to minimize the probability that the 
atmospheric statistics would change between runs. It was hoped that 
in this way the effects of aperture averaging could be separated from 
the effects of nonstationary atmospheric statistics. 

Preliminary analysis c£ the data has indicated definite aperture 
averaging effects. Thus far, however, insufficient data has been 
analyzed to fully verify the accuracy of the theoretical prediction 
of the aperture averaging effect found in the literature. 

We have written computer programs with which to calculate the 
aperture averaging and to compute the log-amplitude variance allowing 
for these effects. These calculations will be completed in the near 
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future . 

2) Probability Distribution. Both Fried [18] and Patton [19] have 
reported departures from the log-normal distribution at 10.6 microns. 
Fried has attributed these results to instrumental errors in the measur- 
ing equipment. It has been, therefore, one of our primary objectives 
to either verify the log-normal distribution for the scintillation of 
a CO^ laser beam or to offer definite evidence that the fluctuations 
are not log normal. 

Our preliminary results seem to favor a log-normal distribution 
but unfortunately they are not conclusive. Out of the first 31 runs 
analyzed, 22 fit a log-normal distribution, 5 a normal distribution, 
and 4 fit neither. Much of the data analyzed has chi square values 
for a log-normal distribution that are well within the 95% confidence 
limits customarily used as a criterion for a fit. That is, the prob- 
ability that random data would have fit a log-normal distribution as well 

is less than 5%. The cumulative probability distribution for two 
such runs is shown in Figures 6 and 7. Figure 6 clearly shows that 
the data fits a log-normal distribution better than a normal distri- 
bution. In Figure 7 the distinction is not as clear. In fact from 
the curves one would have difficulty in deciding between a normal and 
a log-normal model. The value of chi square for the log-normal model 
was 65 as compared to a chi square of 256 for the normal model, clearly 
indicating a better fit to the log-normal curve. In Figure 6, where 
the distinction is much greater, the chi square values were 90 and 
1027 for the log-normal and normal models respectively. 

Some of the data, such as shown in Figure 8, shows a much better 
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fit to a normal distribution than to a log-normal one. Although 
farm conclusions cannot be drawn from the limited amount of data so 
far investigated, it seems that this occurs when a larger aperture 
is used; A possible explanation for a normal distribution could be 
that if the aperture is larger than the correlation distance then 
the intensity incident upon the detector sum of uncorrelated intensity 
fluctuations across the aperture. Involving the central limit 
theorem we may argue that if the aperture larger than the correlation 
distance the distribution function should be normal while for 
apertures small compared to the correlation distance it will be log- 
normal. For intermediate size aperture one would expect a distri- 
bution somewhere between normal and log-normal. It is hoped that the 
complete analysis of all the available data will clarify this point. 

It has been found that an occasional run will have a distribution 
which departs radically from both the normal and the log-normal 
model. Often these runs are characterized by a large negative 
skewness. No explanation can be offered at this time. 

3) Structure Constant . The atmospheric structure constant has 
been computed for a number of the runs, neglecting the finite aperture 
size. The values obtained are within the range of values expected 
from theory. These values will not be reported at this time since it 
is felt that they will be effected significantly by the aperture 
averaging effects. 

Continuation of the Project 

As has been emphasized above, the results which we are reporting 
at this time are based upon analysis of only a small part of the 
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available data. A proposal has been submitted for extension of this 
project. It is expected that when the analysis is completed the 
available data will yield (1) verification of the distribution function, 
(2) accurate values for the atmospheric structure constant, (3) power 
spectral densities, (4) confirmation of the theoretical predictions 
of aperture averaging and (5) information concerning the variation of 
the structure constant and spectral density with time of day and with 
weather conditions. 

The preliminary analysis has indicated that there is every reason 
to expect that much useful information may be extracted from the 
available data. 
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FIGURES 6-9 

The cumulative probability distribution for the illuminance and 
the log amplitudes showing (Fig. 6) a typical small aperture run in 
which the scintillation is log-normal, (Fig. 7) a run in which the 
scintillation is intermediate between normal and log-normal, and (Fig. 8) 
a run in which it tends toward a normal distribution. Figure 9 shows 
an example of an occasional run in which the data departs radically 
from both the normal and log-normal distributions. 
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(3) 

NCI - 

(4) 

NOSCAN- 

(5) 

N0CHAN- 

(6) 

TOLI - 

(7) 

T0L2 - 

(8) 

T0L3 - 

(9) 

LI 

(10) 

L2 


Specifies the number of class intervals to be used in 
the statistics routine. 

Number of multiplexed variables contained in the record 
that belong to a given data set. (For example, 5 multi- 
plexed variables, 200 points per data set). 

Number of variables multiplexed 

Sets tolerance on base limits for base point selection. 

Sets tolerance on signal limits for signal point selection. 

Number of signal base, and amplitude points required for 
a cycle to be completed. 

Number of sequential searches for points allowed before 
cycle is aborted. 

Number of base points required for projected signal point 
search. 
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APPENDIX A 


This appendix further describes the computer program employed to 
reduce the atmospheric data. Included is a list of the FORTRAN state- 
ments of the program. 

The program reads data from digital magnetic tapes with the use 
of three special FORTRAN subroutines, NTRAN, MOVE, TRNSL. These sub- 
routines are used together to convert the data format on the tape to 
a workable format for FORTRAN processing. 

The data is written on the tape in blocks called files. Each 
file contains 300 records of data plus an identification record at 
the first of each file. This record consists of 12 tape characters 
in a binary called decimal format. The 12 characters contain the 
following information; tape number in characters 1-4; scale factor 
in 5; type file in 6(0 = data file, 1 = calibration file); number 
of variables multiplexed in 7-8; number of scans in 9-12. The 
first 4 characters in each data record gives the time. The re- 
maining 2000 binary characters in the record yield 1000 numbers in 
a multiplexed format, i.e., every fifth number belongs to the same 
data group. The program unpacks and stores the desired group in 
an array to be used in processing. The first file of each tape is 
a calibration file. 

The program requires 10 control parameters to be read in as 
data at the first of each tape run. These are: 


(1) IBTYP - Determines method to be used to calculate the signal 

points. IBTYP = 0 will have the program calculate 
amplitude points by taking the difference between the 
average of the base points of the cycle and the signal 
points. IBTYP = 2 will have the program to calculate 
2 amplitude points by subtracting the first signal 
point from the base point preceding it and the last 
signal point from the base point following it. 

IBTYP = 3 will cause only 1 data point to be calcu- 
lated by subtracting the first signal point from the 
base point preceding it. 

(2) IFLAG - The value of this parameter determines the type chi 

square test to be used on the data. If IFLAG = 1, no 
test will be performed. If IFLAG = 2, the data will 
be compared to a normal distribution. IFLAG = 3 will 
call for a log normal test, and IFLAG = 4 will call for 
both log normal and normal tests to be performed. 



INTEGER FILES, FILCk 

COMMON P(IOOO) , AUXI 1000 ) , INDOI 10 ) ,NCNT 1 300) , JOVF, JUNF , 

»HA6H2, 101 ,5im 10) . im. ~ U Ci 2.101 , KK,HH»XL'90,XCID 

COMMON/ AAA/ I PRINT 

C0NHUN/RLL/1RIM 

COMMON/RAT/IRCENU 

CfflWW/RR/ IE(7.30in 

COMMON/UUU/NPNCH 

DIMENS I ON TBUMlW<n 

DIMENSION 1A(3),IB(5) 

CALL KIUm.IJ.l'BI 

IF!I1>155, 150,155 

1S3 PRIN T LS F .ma 

157 FORMAT!* **• END OF FILE •***,181 
“TWIT! 121158, 161, 158 

158 PRINT 159 

T55 — hUR H A I C — *** — READ ERROR — WF*! 

161 CONTINUE 

NTAPt-IBm 

PRINT 160,18 

16U FUKHA T T ' T A P E ~ NUHBbtl-MAA, 

1* IDCAL** ,1AA, 

P UASt-'.lALi 

1* NUMBER OF CHANNELS** , IAA, 


3350 FORMAT I 2110) 

C "READ INSTRUCTIONS AND' TOLERATES 

READ 15,360) 18 TVP, IFLAG,NC 1, N3 SCAN, NOCHAN. LI, L2 

360 F(jkhAT(71 10) 

PEA0(5* 1991) TOLL , T0L2.TDL3 

lyyi” FURMAT ( JE10.0) 

PRINT 361, iarrP,lFLAC,NCl, NOSCAN, NOCHAN, TPLl,T012,T0i3,Ll,L2 


****** ***** 1BTVF 


****** (FLAG 


****** NOSCAN 
"P — NUMBER OF 
****** nochan 


*••«•» toli 

TP rOTETnOK 

****** J0L2 




****** T0L3 


****** l_l 


• irl’f w .1 -I.M :¥ V 


I -j 


M -T l*| i 


****** L2 


511 REAO(5,882)FILES,1START,1CHNO, 1RCEN0 


FI 0,3/, 




1 |J B'l ^ iM J|T; 


k* rM'nMji] | ’ i ht -» t tj t | 



FORTRAN IV 

G LEVEL 

1. MOO 4 MAIN DATE - 69290 13/17/51 

PAGE 0002 

0031 


DO 8892 IJ-1.7 




UO 8892 31*1,300 


00 33 

6892 

1EIIJ.JI 1-0 


CCT3* 


IRECN1 -Isf ART-l 


0035 


IEXIT-0 


5536 

557“ 

1REC-0 


0037 

882 

F0RNATI413) 


SoJS 


PRINT 901, FILES, ICHNO, 1 STAR! , I RCENO 


00 39 

901 

FORMAT 1 1H1* * PROCESSING FILE •**••••••••,15/. 


*• CHANNEL ****•»****••©*♦♦*♦•*••*******•, 15/, 

• • BEGIN AT RECORD ♦*•♦•••***♦*•♦*♦•****•, 15/. 



• ♦ STOP PROCESSING AT RECORD ***»***♦*♦♦• ,15) 


0040 


IRTN-1 


55*1 


1FIF1LES)991,10,10 


004? 

10 

IFIFILES.UT. FILCH) GO TO 9 


oo*5 


CALL NTRANll.lO) 


0044 


GO TU 9 


0045 

991 PRINT 992 


0046 

992 



55T7 


Stoi> 


0048 

9 

CALL RLOAECI F 1LES, 1REC.FILCK.1 BUF, TIME , N* NOSCAN, NOCHAN, ICHNO! 



C 

MOVE TA^E TO OtSlRED FILE 


0049 


IF I F I LCK-F 1LES ) 9,507,097 



c 

READ RECORD, SfORE RECORD IN ARRAY IBUF 


0050 

a 8/ 

CALL REOREC I F IL E S, IREC. FILCK. IBUF, TIKE, N, NOSCAN, NOCHAN, ICHNO) 



c 

MOVE TAf>E TO DESIRED RECORD 


0051 


1FI IREC-1START>8B7,«S9,889 



c 

STORt CONTENTS _ OF IBUF INTO MIN ARRAY P 


0052 

889 

00 88R L»l,N 


0053 

888 

P1L>-1BUF<LI 



C 

READ NEXT RECORD INTO IBUF 


5354 

c 

CALL RED R£CIFilE$,IAEC,FILCK, IBUF, TINE, N, NOSCAN, NOCMAN, ICHNO) 
STORE CONTENTS OF IBUF INTO AUXILIARY ARRAY AUX 


0055 


DU 500 L-1,N 


0056 

500 

AUX1U-IBUFIL) 


005/ 


GO TU 501 


0056 

5u2 

DO 503 L-i,N 


5559 

503 

P<L)»AUX(U 


0060 


CALL RED KECIFILES.IREC.FILCK, IBUF, TIME, N, NOSCAN, NOCHAN, ICHNO) 



C 

STORE IBUF INTO AUX 


0061 


DO 504 L-L.N 


0062 

5u4 

AUXIU-IBUM L» 


006 3 

501 

CALL AVG(N,Ll,L2.IRECNT.TULl.r0L2.T0L3.NCI.lBTYP! 



c 

STUP PROCESSING ON DESIRED RECORD 


0064 


IF ( 1 G X I T > 661*661,507 


0065 

661 

J F 1 JHEC-lftCEND) 502,221,221 


0066 

221 

DU 2293 l«l*N 


0067 

2293 

pm*Auxin 


0066 


1CX1T-1 


0069 

C 

GO TG 501 

PRINT CkAOK TABLfc 


5575 

60* 

CALL PRINT 11! 


00/1 


CALL STATINTAPE, ICHNO, FILES, NC I, IFLAG) 


0072 


GO TU 511 


0073 


FND 





njmmH'TV 

G TE VET 

■ IV'HUir? STB DATE - 69295 

13/17/51 

PACE Ooo 1 

0081 


SUBROUTINE SV'fflNI<,LL,L2,lR£C,r0Ll,TaL2,T0L3,f*CI»ielYPI 



0002 


COMMON Pt 10001 ,*UX( 1000), 1NDQ1 10) ,NCNT ( 300) ,JOVF, JUNF. 





*B'A5El2,101,STG710),IbLOCT2,10) ,l,N7*L90,XLTd 



0003 


COMMON /AAA/ 1 PRINT 



ooo? 


TWHON/RLC/IRTN 



0005 


CDMMDN/PR/IEI 7,300) 



0005 


■TOHHOW/CBA/IBCNT! 21 , 1SCT 



0007 


I REC* I R£C*1 



onus 


151*15) ♦ 1 



0009 


INP*0 



COTTJ 


INDX*1 



oon 


GO TO 1 150,11. 12, 150, 7761.IRTM 



00T2 

130“ 

I RUN** 1 



0013 


IRP-0 



“ours 


isr-o 



0015 


L-l 



rou 


“M»20 



0017 


IF! IRTN.E0.4) GO TO 209 



TOTH 


l)U 1U5U 1-U'JOO 



0019 

1050 

NCNTIIi-0 



— — -UV7U 


j avF-a 



0021 


JUNF-0 



UU72 


“NUAI A*0 



0023 

209 

CALL LIMIT! TOLl.TOLZ, WO 




— C 

SEARCH FUR "FIRST BASE POINT- -IOLB VALuE OF INBX 



0024 

776 

00 1 1»L.M 



(JU75 


IFII.GT.NK! OUTD 2000 



0026 


IF!PH)-XL90)1»3,3 



0527 

3 

INDX*1 



0020 


GO TO 4 



8023 

j 

enMTIKUE 



0030 


IE11, lR6C)*ieil,lRbC)4l 



oon 

c 

‘IRTW*4 

RETURN TO DRIVER FOR NEXT RECJRD 



0032 


RETURN 



0033 

2000 

L* L 



003? 


'H-H-TJR 



0035 


IRTN-5 



— — 0035 


'RETURN 



0037 

4 

M-1NDX415 



0030 


"L«IHDX 



0039 


CALL LIN1UT0L1,T0L2,NK) 



oonj 

/ 

TSC-I 



0041 


IBCNTCD-0 



- otm 


"IBCNTI21-0 



0043 


ISCT-0 



50?? 

■" 10 

LC*0 



0045 


NP-0 




— c 

BASE POINT SEARCH 



0046 

11 

1 F ! P! I NDX )-XL90) 12,13,13 



00?7 

“13 

I BCNT ( 1 BC J* I BCNT t 1 BC ) 



0048 


1 * t BCNT (IRC) 




C 

STORE BASE POINT 



0049 


BASE! IBC,I)*P1 INDX) 



8038 


■iBincriK.n-iNDX 



0051 


IF11BCNTI IBO-IO) 14,401,401 





— 



1, MOD 4 AVG OATC * $92?Q 

13/1T/St 

PA6E 0002 

« 

c 

— 






• 

— 

5351 


IEt2,IREC)»IEC2,lRECIH 



Ml 

( 

3555 

Pi — 

1N0X-INDX ♦ l 



• 

0 — 

5357 

158 

IRTN*2 

NORMAL RETURN TO DRIVER FDR NEXT RECURD 



II 


— sws 


RETURN 

SIGNAL POINT SEARCH 




\ > 

3359 

12 

1F(P( INDX1-XL101 20,20,15 
IF1LCI 16, 16,21 



1 n 

0 ~ 

0641 

16 

NP-NP*1 

CHECK NUMBER OP UNSUCESSFUL PASSES 



'2 


0062 


lf(NP-U) 1944 • 1944, 20, 
1 F 1 ISC7.E0.0) GO TO 14 



i 

V 

0044 ' 


zrmrTo 

r F t P 1 1 NQX I-XL90) 22,30,30 




0 — 

0066 

30 

IFtlRUN>l0,10,31 

SFT BASE COUNT 1NDX CONDITION 





0067 

31 

IF t IBC-i > 34,33,34 
1 8C* 2 



> 

U 

— 0060 

34 

GO TU 10 
I BC* 1 




0 "" 

— 0671 " 

20 

c 

1 F ( I RUN 1 21,23,24 
CHECK LOOP CONDITION 





0072 

23 

IF ILC) 60,40,29 




u 

5377: 

25 

lOl 




0 “ 

0076 

29 

c 

1 sex » i set ♦ 1 

STORE SIGNAL POINT 



) 


5577 


SlfitlSttl-PlINDXI 
1N0Q( 1SCT )«Ii4DX 



> 

1 

0079 


IF ( ISCr-10> 26,501,501 
IRTN-4 




o 

0081 


"ifei3,lR€C»-iei3f IRECJ41 
RE TURN 



) 


008 5 

26 

INDX*1N0X«1 
I F ( INUX-NK I 12,12,162 



> 

V 

0085 

162 

IRTN»3 

NORMAL RETURN TO DRIVER FOR NEXT RECOTO 




0 “ 

0086 

22 

RETURN 
NP*NP+ 1 



) 


ooea 


IF(NP-L1I 26,26,208 

r.Hfr.K NUMBER OF BASE AND SIGNAL POINTS 



0 

u 

0089 

40 

IFI IBCNT 111. LT. TOLU GO TO 222 
(FI IBCNT12 I.LT.T0L31 GO TO 222 




0 “ 

0091 


iFUSCT.LT.TUL3) GO TO 224 
GU TU 100 




13 

0093 

222 

IEI4*IREC)*IE14, 1R£C)*1 
GO TO 555 



‘ 

U« 

* 

0095 

224 

lE(5.IRECI-Ifc(5,IREC)Al 
GO TU 555 




c’ ~~ 


C 

ENTER ERROR RECYCLE 



) 

B 

0090 

555 

L * INOX 



0 

\ '4 — 




run IRAN 

iv o Lc vc L L, MOD 4 AVG DATE * 6$290 

1J/17/51 

PACE 0003 

0099 

0100 


M-IN0X420 

IRUN-1 



0102 


1M 1PK1NI GO TO 742 

PRINT 666, IRECUNDX 



0104 

666 

FORMAT ( fX» fttC" ' , I 10 » 10X, * INOX** 1 1 1 0 1 
IFUREC.eO.IRP> GO TO 742 



0106 

333 

PRINT 333* I P C 1 1 » I * 1 f NK I 
FORMATUX,/, tlX,lOF10.3> I 



0108 

742 

IKK-lKtO 

CONTINUE 



0110 


IF C INP~5 >209* 921* 921 



UI1I 

0112 

' — 92 1 

”11:1 / • IKtO-IEI 7| IHtC)+l 
IRTN - 4 




c 

POINT SEARCH CYCLE COMPLETE 



0115 


~~iuni*u.o 

SUM2*0.0 




c 

PREPR0SSIN6 FOR HISTOGRAM FOLLOWS 



0118 

101 

UU IU1 J”i,I — 

SUN1-SUM1* BASE 1 1* J 1 



0120 


00 102 J-1»I 



0122 

IU2 

5UM2*SUM2 +6AStf^*J) 

BA* ( SUM1+SUM2) / 1 1 BCNT (11 ♦IBC<T(2)> 




c 

LALL HIST (¥A7Nr , NC I P IBTYP* I RUST 
SET BASE COUNT INOX C0N0IT10N 



0125 

43 

45 

ITU BC-1J 44*45 * 44 
1 8C*2 



UI^6 

0127 

44 

SO TO 46 
I BC*1 



0129 

46 

1 5CT*0 — " 

IBCNTUBO-O 



0131 


LC* 1 — — 

NP«0 



0133 


IRUn-o 

IN0X2-IN0X*. 



0135 


L-INOX “ — — 

H-1N0X2 



0137 


U*LL LIMiV T IUL1, TDL2*NK) 
GO TO 12 



0138 

70 

LUO* INTO NEXT RfJttikb FOR S 161 At POINT 

IF|1BCNT(IBC)-L2) 14,72.72 



0140 

72 

"C«Tndx — 

M-INOX ¥ 9 



0142 


LALL L 1 M I T 1 I UL 1 * TUL 2 * NK J 
IF(PdNOXl-XLlO) 20,20,14 



tP*U — — 









FORTRAN 

IV 0 LEVEL 

1. MOO 4 LIMIT DATE * 69290 13/17/51 

PACE 0001 

1 

0001 


- C 

SU6K0UT1NE LIMlT(T0Li*T0L2«NK) 



5535 


COMMON PI 1000 ) , AUK ( 1000) , 1ND0I 10) ,NCNT 1 3001 , JOVF, JUMP , 
♦BASE 12. 10 1 ,51010), 1 BLOC (2, 101 . KK.MM.XL90. XL*10 



0003 


1K-KK 



0004 


IM-MH 



0665 “ 


10*0 



0006 


1F1IK-NK) 502,501,666 



0007 

666 

POINT 3000 



0008 

3000 

FORMAT (• IK GREATER THAN NK •) 



0009 

501 

AMAX-Pl I K 1 



0010 


AMIN-AHAX 



cron 


I 0» 1 M-NK “ 



0012 


GO TO 361 



0013 

502 

1EUM-NK) 360.360.361 



0014 

361 

10*1 M-NK 



55T5 


IK* NX 



0016 

360 

AMAX-PUK) 



55TT 


AHlN-AMAX 



0018 


DO 350 



0019 


IFiAHAX-PI J) > 301,302,302 



0020 

301 

AHAX-PUJ 



5021 

302 

IFIAHlN-PljT) 350,303,303 



0022 

303 

AMIN-P(J) 



~TO77 

3*0 

CONTINUE 



0024 


1 F t 103 380,380.381 



5557 

381 

AMAXP-AUXl 11 



0026 


AH1NP-AMAXP 



5557 


DO 650 J-1,10 



0028 


IFUHAXP-AUXUH 601,602,602 



0029 

601 

amaxp-auxTj 1 



0030 

602 

1F< AHINP-AUXI J) ) 650.603,603 



Soil 

603 

ANINP-AUXI J) 



0032 

650 

CONTINUE 



5573 


1F(AMAX-AMAXP(700,701, TOl 



0034 

700 

AMAX»AMAXP 



0035 

701 

1 F < AM IN- AM IMP >380,360) 703 



0036 

703 

AMIN*AM|NP 



5377 

380 

A- AM AX- AMI M 



0038 


XL90»AMAX-T0Ll*A 



' 00J9 


XElO-ANlN* T0L2*A 



0040 


RETURN 




66*1 End" 



f OUTRAN IV V LEVEL 1. 


HOD 4 


MAIN 




DATS 




0002 

0003 

oS5i 

0003 

TOST 

OQOT 

TOST 

0001 

0010 

TOTT 

0012 

TOT 

0014 

TOT 

0016 


TOT 

0019 

tstot 

0021 


0022 


0024 
002 9 
TO* 

002J 


0028 

0029 

003V 

TOT 


0012 


0013 


0034 

0035 
003* 
0031 

0036 

0039 

0040 

0041 
0042“ 


0 


c 


z 


c 


T 


c 


c 


statistics 

SUBROUTINE ST AT (NT APE tNCH iNFIL E *NCI * IFLAG) 

J El AG « I NO CHI SQUARE TEST 

2 Chi square test on normal distributi 

3 CHI SQUARE TEST ON LOO NOftHLL DISTRIBUTION 

4 CHI SQUARE TEST ON BOTH 

C O HHUH 1* 1 1 CPU ) , AUX 1 2 OOP! . 1 NODI 1 0 > .HCNT1 1 300 ) , JOE ■ JUf , B AS6 1 2 . 10 1 . 
4PIG1 10) , I810CI 2.10 ) ,Rk7rM,XL93,XL10 

COMMON/ UUU/NPNCH . 

DIMENSION YI300),XLNI300),Q<4> ,RI4) 

DIMENSION NCNTOOOI 

OO 60 1-1.300 

60 NCNTII )-NCNTl(I) 

C-1000/NCI 

F I ML) THE HIGHEST AND LOWEST CLASS 1 NT CM VAC 

00 1 i-wnCi " 

1 FINCNT 111) 1,1,2 

2 1L0»1 

GO TO 3 

1 CONTINUE 

i yu 4 |. ice , nci 

If INCNT I J ) )4,4,5 

5 1H1-I 

~ 4 CUnTiNUE 

FIND NUMBER OF POINTS AND FLOAT NCNT 

,1«0 

00 6 I-IL0.IH1 

Ylll-NCNIin 

6 N-N4NCNT 111 

PRINT HtAOINCS I 

WKIT6I6.101) NTAPEiNCHtNFILE 

MRlTLlbi 102) 

Defalt due to too few class intervals 

NN» 1H1-1L0 

I F INN- 10 1 50.50.51 

SO nR1TEI6,110> 

return 

FIND AVERAGE AMPLITUDE 

91 XU»fi 

AVE«(T.OO 

DU U I-1L0.1HI ___ ' 

Xl»I 

0 AYE-AVb»V( 1)41X1-0. SI *C 

COMPUTE CUMULATIVE PRUBADI1 IT l EC AND LOG AMPLITUDES 

SUM-0.00 

DU 10 l-ILO.IHI 

XI-I * 

xi»ixi-o.s)*c 

XLNI I >«0.5*ALOOIX1/AVE ) 

SUM-5UH4Y 1 I > 

CP-SUh/XN 

ID WIUlh<6.I03) Xl.XLN Ul.Ylll.CF 

M>UTEI6tlll) JOF, JUF 

COMPUTE HQMLNTS ABOU T THE MEAN 

XLA-0.00 


■-J 

t— 1 


0043 



Q!4>-Q(J)*Y( l)«UXI-Q.5)*C-AVE >»»J 


^aiwinimrviiinxiwaiiu Til 


u a A3 J>2 f * 


r*immnmn 


2) RIJ)«RIJ)/XN 


IDIIIQ iii[] 


WR1TE<6, 10*) NT 


M I *T‘ > | (U fj I ■ 


SIGL-SURKRI2) ) 




SKL - 0.5»R(3>miSL**3) 


M tii'rrtuf iruttMt t yUPWtlfi’f 


XXURL » (IR<4>/<RI2)*«2)l-3. 01/2.0 


rrwn. nuncnu 

WRITE (6, 105) AVE,SIG,SKEW,XKU1 ,XLA,SIGL,SKL.XKURL.N 


■ ()■,) | ,1 


PUNCH WOO. NTAPEtNCHtNFILEf AVEt SIG* SI GL t XL A 



chi square tfst ~ 

32 CALI CHIICSQ.Y.ILU.IHI.C.NUSE.AVE.XLA.SIG ,XN,0.SIG) 


Mnl ItrlWfluol L1U|NU>C 

PRINT CHI SUUAAE NORMAL 


JU IU UlfJ • 

42 CALL^CHI (C50,Y,1LQ. IHI .C.NUSt. AVE ,XLA,SIG ijjljilitltL) 

PRINT CHI SQUARE LOC NOMAL 


KC IUHN 

10) FURMATl'l'.SX.'IAPE NUMBER '» A4. 5X» 'TRACK' t I3.5X. 'FILE' ,13) 


| ( 1 7 *^ Jlj.Kf.ll 


I ri WXTltJw |iT tjd 


JB | £f '>] B 


1 BX. 'CUMULATIVE PROBABILITY'/ I 


UmilLYI 


104 FORMAT (' O' v IGX« 1 NUMBER OF CLASS INTERVALS «',16) 


JinT.T Vi^:y\H 4 Jtr> ftlif ' 17 VI ibi] d’l fig inn 1 PHn 3 T !Mi'tiB 


*13X,'KURrOSJS'/7X > 4E21.6/7X,4E21.6/'0' ,10X,'NUMBER OF DATA POINTS' 


106 FORMAT! 'O'iIOX.' CHI S0UARE-' ,E 14. 6/ 1 IX • NUMBER OF CLASS INTERVALS 

lUjtd _ I “1 5| 

110 FORMAT I >0' *5Xi 'TOO FEN CLASS INTERVALS •/ 

r-H>r;5Xr — ETFCuT fM _ 3J^S TATTST ICTTAl CO lTTTON S0fmi0£jTJ 

111 FnRMATI *0* ,5X, 'NUMBER OF OVERF LOWS' , 16/ 

*6X, 'NUMBER Of und£RFlom&' ■ 15 I 
FN0 






COOL 

0002 


UTOT 

0004 

oob5 

0006 

SWT 

0008 


oooi 

"TOTo 

0011 

0012 

0013 

0014 


TO IT 
0018 
fforr 
oozo 
obiT - 
0022 
0023 


"00 Z 5 
owt 


002B 

0021 

obSb 

0031 

1032 

0033 

0034 

0036 

75516 - 

0037 
003d 
0039 


CHI SUUARC TEST 

SUBROuf 1NE CHt ICSQ,Vl,ILO,IHl,C,NUSE,AVE,)(LA,SD,XN,NTYP,SX> 

r> i he ns ion yi3oo.3) , vnaoo) 

UU 1 I«l,300 

1 YH.1L0.00 

.MUSE-0 

KH«1 11U41HII/2 

J-l 

YU,2)-»VC-10.Q»S0 

"E GRUOP CLASS INTERVALS ON UM END 

00 2 1- ILO.KW 

yij.d-yu i ) ♦ riJii) 

lFIVU.n-5.01 Z,Z, 3 

3 YU.31-CM 

NUS6 » NOSE * 1 

J - J ♦ l 

vu,. 2 ) - c*i 

2 CONI HUE 

_C GROUP CLASS INTERVALS ON HI 63 SIDE 

, , , H j 

YIJ,3t - AVE ♦ 10. 0 «SU 

6 IeU-KHI 10,10.4 

4 Y1J.13 - YU, 11 ♦ Y1U1 

IFt»lJ,l)-S.(H 11, If. 5 

6 Y(J,21 ■ C«ll-ll 

'MUSE - OUSE ♦ 1 
J - J » 1 

yu, 3) • fSTFir ' ~ 

n i * i-i 

OU fO 6 

± COMPUTE THEUHIT1CAL PROBABILITY 

10 CbU • 0.00 

UU 30 1-1. MUSE 
XLL , y(i,2) 

KUL.YH.3I 

24 CALL STNPT FTH.XLLtXUL./l.NTYP.AVE.SX.XLAJ 

C COMPUTE CHI SQUARE 

1F1ETH) 31,31,30 

31 HR 1 TE 1 6, 1001 I 

100 FORMAT ( * 0 ZERO VALUE UF THE0R1 1 1CAL PROBABILITY IN', 16,' TM 
_ 1 INTERVAL'/ 6K, 'EXECUTION OF CHI SQUARE TEST DISCONTINUED' 1 

RETURN “ 

10 CSu»CSU'IIY(l, 1)-XN»FTH)*«2I/T XN4FTH1 

RE TURN 

EN1) 



0001 



SUBROUT lilt SIHP<SUH,FLL.FUL.N,NTVP.A,».CI 


r 


iNtEGRANO FUNCTION REMOVE XHENCHANGING FUNCTION 

0002 



PBFtX.A.SI » U.O/IS*SORTI6.2a3ia>))*EXPI-O.S*(IIX-A>/S)**21l 




FNP.N-1 

0004 



L)6LX*» ( FUL-FLL > /FNP 




SUM-0. 0 

0006 



SUMl-u*Q 

5OT7 



SUM2-0.0 

0008 



0U l 1-ltN 




FK-t-1 

0010 



X*FK*DELX+FLL 


— r 


CALL FUR INTfcRGRAND SUBROUTINE HERE 

oou 



IFINTVP) 101,101,110 



101 

VAL*PHFU*A7ei 

0013 



GU TU 102 

CBHi 


"TTo" 

IFIX) 20,20,100 

0013 


20 

VAL-0.00 

0016 



GO ru 102 

001? 


100 

XX*0.t>*ALOGU/A) 

fflMl 



V XL *PPF 1 XX , C ,B ) 

0019 



VAL«0.S»VAL/X 

0020 

c 

102 

CONI INUE 

HfilBliItM ■■ 



1F< Utg.l.UR.l ,£U.N) GO 10 2 

0022 



J-PUDI1.2) 




IF(JI3,*,3 

0024 


3 

SUH1-SUH1+VAL 




GU TU 1 

0026 


4 

SUM2»SUH2*VAL 

00 2 ? 



GO 10 1 

002b 


2 

SUH»SUM*VAL 

0029 


L 

cuNUNue 

0030 



SUM-SUM+2 • 0*SUH1 ♦ 4 # 0*SUM2 

00)1 



S UH- SUM* OL LX/3.0 

0032 



RETUKN 

0033 



tNU 


0002 


SUBROUTINE KlOllJ.U.ICi 

1 NfEGER 16tS),IA<3), BUFIS0D.FLCNT.TTB13I .BLK ,f 1 LES ,F 1 LCK 



14 ICUI-1BU) 


Nad PM«n _ 


CALL NTMNIL,22> 


10 TIKE'BUFIl) 




*»3+2»lCHNO 


IBUEtD-O 


CALL HQVc ( t BUE 1 1 ) >K» 1 1 
















FURTRAN IV b LEVFL l, MOD 5 


t*R 1 o T 


DATl- » 69290 


13/17/51 


PAGE 0001 


OQOi SUBROUTINE PRINT I NNN ) 


000/ 

0003 


CUMMUN/PR/ I £ ! 7 , 300 ) 
COMMON/RAT/IRCEND 

5UC? 

0005 


DIMENSION IESUHI7) 
PRINT 3 


3 

F0KMAT(lHl t 50X,'0ATA PROCESSING I RREGUL AR T IES • , / , 10X , ' ERROR CODES 
• FOLLOW,////, 


* 10X f f NO BASE POINTS FOUND IN 3 ASE SEARCH ****•*•*•• 1'/, 

* IPX , * NUMBER UF BASE POINTS EXCEED S 1<L »♦♦»»»» »»»♦♦ 2'/, 

* 1 OX , ' NUMBER UF SIGNAL POINTS EXCEEDS 10 *********** 3 »/, 

♦ IPX, 'NUMBER UF BASE POINTS INSUFFICIENT ********** » A'/. 



•10X, 'NUMBER UF SIGNAL POINTS INSUFFICIENT ***•**••* 5V, 

•10X, 'NUMBER UF PASSES EXCEEDS LIMIT LI •**•••*••••• 6'/. 



•10X, 'NUMBER UF ERRORS IN RECORD EXCEEDS 5 ********* 7») 


0007 

UU 555 1=1,7 


oooS 

555 1 ESUMI 1 ) =0 


0009 

PRINT 1 


trcTo 

1 FORMAT! 10X,' ERROR', HX , * l ' , 1 3< , '2 * , 13X , ' 3 ' , 1 3X, '5 ' ,13X , • 5 • , 13X, 

♦ '6' , 13X, *7' ) 


OoU 

PRINt 100 


001/ 

100 FURMATUOX.'RECURD',///! 


0013 

DO 65 l * 1 , I RCEND 


0015 

DU 55 K=l,7 


0^13 

65 1 ESUM! K 1*1 ESUMIRl + IE t K , 1 ) 


0016 

DO 2 1=1,1 RCEND 


0017 

00 12 K= 1 » 7 


0010 

IFIIEIK.D) 15,12,15 


0019 

FS continue 


0020 

r.U TO 2 


ooIT 

13 RR1TE16.5) 1,(]E(K,1),K=1,7) 


0022 

5 FORMAT 1 10X, I 3, 10X,I5,10X,I5, 13X,I5,10X,15,10X.I5,10X,t5,10X,I5l 


0023 

2 CONTINUE 


0025 

PRINT 20 


0023 

20 FORMAT! lux,//, 10X, 'ERROR CODE' , 10X, 'NUMBER OF ERRORS') 


0026 

DU 50 K*l, 7 


002 7 

50 PRINT21.K, itSUM(K) 


002 B 

21 FORMAT t 10X , 16, 15X , 19) 


0029 

RETURN 



0030 END 
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APPENDIX B 


This program generates a set of N random numbers having a log- 
normal distribution and a pre-selected mean and standard deviation. 
The program is in the form of a FORTRAN IV subroutine. 

Theory : By definition a log-normal random deviate is one 

whose logarithms are normal random deviates. Thus if (X.) is a set 
of log-normal random numbers then their must exist a set of normal 
random numbers (y ) related to the X^ by 

y. = In X. Bl 

i i 

Equation Bl may be generalized by the addition of appropriate scaling 
factors, i.e. we may let 


y. = a In X. + b B2 

i 


Now by choosing the mean and variance of the (y.) and the values of 
the scale factors a and b it is possible to generate a set of (X.) 
having any desired mean and variance from a set of normal deviates 
(y^) . Solving B2 for X^ we have 

y± ~ b 

X. = exp ( ) B3 

l r a 

Since we wish to specify only two parameters, viz., the mean and 

standard deviation of the (X.) it seems reasonable to assume that 

we will need only two parameters in equation B3. We therefore let 

a = 1 and take the mean of the (y.) to be zero. B3 then becomes 

■'l 


= exp (-b) exp (y^) B4 

taking the average of both sides of equation B4 we have 


X = exp (-b) exp(y 1 ) 


B5 
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and also taking the second moment of (X^) about zero 

X 2 = exp (-b) exp (2y_^) B6 

the averages of the 
be evaluated easily 


-expCnyJ 


exponential functions in equation B5 and B6 can 


- 1/2 


= (2irt ) 


exp(ny) • exp(y /2o)dy B7 


-x 


Combining equation B5,B6 and B7 we obtain expressions which may be 
solved for the scale factor b and the required standard deviation 
of the (y j .) 


a 2 = In (y/X 2 ) 

and 


B8 


exp(-b) = Xexp(-a /2) 


B9 


where y is the second moment of the (X_^) about zero. 

Program: The log-normal generator makes use of the normal 

random number generator included in the IBM Scientific Subroutine 
Package for the 360 computer. This routine (GAUSS) generates normal 
random numbers with any required mean and standard deviation. Coding 
for the program is shown in the accompanying listing. The argument 
list is as follows: 

AVE The required mean. 

VAR The required standard deviation 

Y A vector of log-normal random numbers returned by the 
subroutine. Y is dimensioned by the calling program 

N The number of random numbers to be generated 
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IX A "seed" required by GAUSS. IX must be a 5 digit odd integer. 

Statements 003 to 006 compute the required standard deviation for 
the Gaussian-random numbers and the proper scaling factor. Statements 
007 to 009 call GAUSS compute a log-normal random number from equation 
B5 . 


The log-normal random number generator has been used to test the 
statistics subroutines used in our data analysis program. A Calling 
Program which will provide the subroutine STAT and CHI with either 
normal or log-normal data is given. 



FORTRAN 

0001 

0002 
0001 
0004 

ooos 

000 & 

000/ 

□008 

0009 

0010 
0U11 


IV G LEVEL 1, MOO 4 LOGN OATE « 69290 14/18/0* CAM 0001 

SUBROUTINE LQGNI AVE , VAR, V.Nil X I 
DIMENSION Yll) 

VAR - VAR**2 ♦ AVE**2 
S1G * AL0GtVAR/AVb**2) 

2BAR-EXPISIG/2.0) 

SlG - SORIISIG) 

00 1 I-l.N 

CALL GAUSS ( IX, SI G, 0,0, X) 

1 Y(ll-IAVfc/ZBAR)*EXPIX) 

RETURN 

END 


00 

l- 1 



FORTRAN IV l» LEVEL !• MOO * 


MAIN 


DATE - 69290 


14/U/09 


FACE 0001 


0001 


COMMON NCNTI300) 

0002 


DIMENSION Yt 5000) 

0003 

20 

READ1 5, 1001 AVE.VAA, N, IX 

0004 


1FUT) 10,11,12 ' 

0005 

10 

STOP 

00G6 

11 

WRITE 16, 1031 

0007 

103 

FORMAT! >1 LOG NORMAL DATA') 

0000 


CALL LOGNiAVE ,VAR*Y,N,IXI 

0009 


GO TO 13 

0010 

12 

WRITE! 6, 104) «* 

oou 

10A 

FORMAT 1 1 1 NORMAL DATA*) 

0012 


DU 14 1*1. N 

0013 


CALL GAUSS I 1X,VAR»AV£*ZI 

0014 

14 

YIII-2 

0015 

13 

DO 2 1*1,100 

0016 

2 

NCNTU) - 0 

0017 


L*0 

0018 


M*0 

0019 


DO 1 I»l, N 

0020 


J*Y1 I >/10.00*l 

0021 


IFU.LT.ll L*L+1 

0022 


IFIO.GT.IOO) M-M+l 

0023 


IF) J.LT.1.0R.J.GT.100) GO 

0024 


NCNTI JI-NCNTI •))♦! 

0025 

1 

CONTINUE 

0026 


CALL ST AT 10, 0, 0,100 *4, L,M) 

002 7 


GO TO 20 

0028 

100 

FORMAT (2F10. S, 3110) 

0029 


END 


c 


00 

t-o 
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APPENDIX C 


The following is a typical computer print out for a file of data. 
The control parameters are listed at the first of each new tape run. 
After a file has been processed, the irregularities found in each re- 
cord are listed in tabular form. The statistical calculations made on 
the data is then printed in a tabular form and identified as to it's 
tape, track, and file number for later reference. 





COOES FOLLOW 


TJO R ASf 
NUMRCR 

: P 
OF 

Ol NTS FOUNO IN BASF SEARCH ♦•******♦♦ 
RASF POINTS bXCb p l>S 10 •'*•***♦**♦*♦ 

1 

2 

n\j natfr' 
NUM3FR 

7TP 

OF 

SIGNAL P(5f NTS EXCEEDS 10 **•**••***• 
BASE POINTS I NSUFFI f I FNT ***♦*#**♦*♦ 

3 

4 

UUIZblidi 

"77 

■4 RCT T |WJi| liMUii M fJTt ^TTTTTT't'TY^^t^ 

NfjMRER 

OF 


6 


QRD LXCEEDS 5 *•*•*••** 
? 



F«f»no chop 


NUMBER OF ERRORS 


fRRFGUL ARTIES 











TAPC NUMBER 3397 TRACK 4 FILE 5 

AMPLITUDE 

LOO AMPLITUDE 

COUNT 

CUMULATIVE PROBABILITY 




0. 11 5000 F 01 

-0.3879976 00 

0.7000006 

01 

0. I91058E— 03 




n.i?bouor "ui 

-0.3463Q6E 00 

0.300000E 

01 

0 • 477646E-01 



' 

i)«nSOOQF 03 

-0*3078266 00 

0.1 200006 

02 

0. 16Z400E-02 




0. 145000? O'J 

-0.7770 5 6E 00 

0.520000E 

62 

0.P59151 E-02 



— 

0.155000E O' 

-0.2387506 00 

0.0800006 

02 

0.1494S1E-01 




0. I45030F' CTC 

-0.207490F 00 

0.1A6000E 

33 

0.308559F-01 




{>.1750006 03 

-0.178070E 00 

0.267000E 

03 

0.543622E-01 




0*Td500GF 03 

— 0*l502fc5F 00 

O.196000E 

03 

0.94I91BE— 01 




0* 1 95000 r 03 

-0* 1 2T96TE 00 

0 .5060006 

03 

0.14Z530E 00 




O.70?OdOF 03 

-0.9895T9E-01 

0.6740 ooe 

6i 

0.2049146 00 




0*21 •jOOOL 03 

-0.751439F-0I 

0.734000E 

03 

0.2770356 00 S 




0.2750flfl£ 03 

-0.524127E-01 

6 .*<60007 

03 

0.3568036 00 




0.735000r 03 

-0.30S702E-01 

0.8490006 

03 

0.4399126 00 




0.2450UfjF 01 

-0 *981 380E-0? 

O.ABloOOE 

33 

5. 5245756 00 




0.255000F 0 3 

0.1016896-ftl 

0.0670006 

03 

O.6OA097F 00 




0* 7 650 IGF 03 

0724AtH8F-01 

0.B0600QE 

65 

0.&83894E 00 




0.775000H 03 

0.479222E-01 

0.6990006 

03 

0.7506A9E 00 




0.7850006 03 

0.6578156-0 1 

OTSAfOOof 

'53 

0.806171E 00 



' 

O.2950OOP 03 

0.830243F-01 

0.5030006 

03 

0.854Z22E 00 




U. lUPO'J’IF 03 

0.9969286-01 

STATfo'ooE" 

<53 

6. *434*56 50 




0.315000= 03 

0.11582 3E 00 

0.2940006 

03 

0.92I570E 00 




0.325000^ 03 

6.I31449F 00 

0.249OOdE 

03 

5.4V3357F 00 




0.315000C 01 

0.1466026 00 

0.I81000E 

03 

0.4620396 00 




0.34 500(16 03 

0.161109E 00 

0.123000E 

03 

0.974589E 00 




0.155 lOD* 01 

0.175596E 00 

0.970000E 

0? 

0.9838566 00 




<1. 1&50UOF 01 

fr.'189486E 00 

O.SrUoooe 

02 

6.9*4 1«4F 00 




n.3 75Qi)0r 03 

0.7030006 00 

0. 460000E 

02 

0.9537916 00 \ 




0.3fl5090 r 03 

0.Z1AI59F 00 

O.i^OoOOF 

02 

0.994370E 00 




0.3O5(| j'jr 03 

0 • 228980E OO 

0.1700006 

02 

0.9979946 00 




0. VO^OOOi 03 

0724U81F 00 

0.100000E 

02 

6.498949? 00 



— 

0.41 bOOOf 03 

0.251A7AE 00 

0.3000006 

01 

0.999236E 00 




0. JOT 03 

0.2655IO.F 00 

0 .3000006 

m 

0.999527E 00 




0.435000r 03 

0.2T721OE 00 

0.500000E 

01 

0. 100 00 0 6 01 





o= fJV co rinws o 

nims r q v rmrnTTO: jr 



NUMBi'i uE CL A 33 INTERVAL? 

Ji- 





AVlR AfJF ST ANOA«l) OtVlAHnN 

0.24 )S>4h 03 0.470818F 02 

Skewness 

0.1 79I44E 00 

KURT0S1S 
0.190 1056-0 1 



-IP. 1976 >7F-J2 

0.950000E-01 

^78171 84F-0I 

-0.471 ifliP-Ot 



NU y l c r r iF HUi PITOTS 

1 0468 





CHI 0.7*irj7l*F 03 

NU«3I K P F CLASS WEPVAI S 

U$Fn * 29 




- 

CHI 02 






NU'IH C F 71 r LA55 1NM.PV4LS 

uafo - n 







00 

a** 
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